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FOREWORD 


PART ! 


This lo o combined report for the two projects at the Massachusetts Institute of 
Technology which arc sponsored by the Office of Naval Research under Contract N5orj60. 


Project on Machine Methods of Computation and Numerical Analysis 

This Project is an outgrowth or the activities of the Institute Committee on Machine 
Methods of Computation, established In November 1950. The purpose of the Project Is (l) to 
Integrate the erforta or all the departments and Jroups aL M.I.T. who are working with modern 
computing machines and their applications, and ( 2 / to train men In the use of these machines 
for computation and numerical analysis. 

People from several departments of the Institute arc taking part In the project. 

In the Appendix will bo found o list of the personnel active in this program. 


Project Whirlwind 

This Project makes use or the facilities of the Digital Computer Laboratory. The 
principal objective or the Project la the application of an electronic digital computer of 
large capacity and very high apecd (Whirlwind 1) to problems In mathematics, science, engi¬ 
neering, simulation, and control. 


The Whirlwind I Computer 

Whirlwind r la or the high-speed electronic digital type. In whloh quantities are 
represented an discrete numbers, arid complex problems are solved by the repeated use of 
fundamental arithmetic and logical (l.o., control or selection) operations. Computations 
are,executed by fractlonsl-mlorosecoud pulses In electronic circuits, of which the princi¬ 
pal ones are (l) the flip-flop, a circuit containing two vacuum tubes bo connected that one 
tube or the other la conducting, but not both! (2) the gate or coincidence circuit-, (5) the 
magnetic-core memory. In which binary digits are stored as one of two directions of magnetic 
flux within fcrro-magnotlo cores. 

Whirlwind I uses numbers or lo binary digits (equivalent to about 5 decimal digits). 
This length wan selected to limit, the machine to a practical size, but It permits the computa¬ 
tion or many simulation problems. Calculations requiring greater number length are handled by 
the use or multiple-length numbers. Rapid-access magnetic-core memory has a capacity of 
32,768 binary dlglta. Present speed of the computer la b0.000 single-address operations per 
second, equivalent t.o about 20,000 multiplications per second. 


Machine Methods of Computation and Numerical Analysis 


1, GENERAL COMMENTS 


During the quarter reported here, Integration of the work of the Project with tnc 
general research of the Institute has proceeded apace. Whirlwind 1 has been used by a wider 
variety of research projects and the Machine Methods and Numerical Analysis staff nave helped 
translate a wider variety of computatIona 1 problems Into machine language than previously. 

One new area of note is In the difficult field of weather forecasting. The Project sponsored 
a series of seminars wnore some of the baala problems were discussed by members or the Meteor¬ 
ology Department and by ;rofeaaor Norbert Wiener. The report entitled "Multiple Prediction 
Theory outlines the resulting work to date 

The research projects reported in the following pages represent a wide variety of 
subjects: their common aspects lie In the adaptation of the problem for machine computation. 
Each application has added to out- experience in 'he techniques of such application: many of 
them hove resulted in useful sun-routlnes, which are now in the Whirlwind library, available 
to other workers. 




» 
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2.1 Index to Reports 


Title 

Multiple Prediction Theory 

Counting Structures or Finite Relstlons 9 

Multiple Scstterlng of Waves from a Spatial Array or Spherical Scatterera 10 

Energy Bands In Graphite 17 

Variational Determination or Atomic Wave Functions 10 

Scattering of Slow Electrons from Atomic Oxygen 19 

Neutron-Deuteron Scattering 20 

Eigenvalues for a Spheroidal Well 21 

A Seir-ConalBtent Determination of the Nuclear Radius 23 

Dynamic Behavior or Shear Wall Testing Machine 2A 

Analysis or a Two-Story Steel Frame Building Under Dynamic Loading 21 

Response of a Five-Story Frame Building to Dynamic Loading 21 

Machine Solution or the Diffusion Equation 25 

Evaluation or the Explicit Difference Formula Tor a Parabolic 

Differential Equation 29 

Stability or Thin, Shallow, Elastic Sheila 31 

Steepest Descent Analysis or Resistivity Data 32 


GRADUATE SCHOOL RESEARCH 


2.2 Progress Reports 
MULTIPLE PREDICTION THEORY 


During the last quarter some work on prediction theory was begun and completed. 

This work waa stimulated by two of the project seminars In which the application, computatlcn. 
and theory or prediction waa presented by Professor Thomas F. Malone and Mr Robert 0 Miller 
or the Meteorology Department and Professor Norbert Wiener of the Mathematics Department. The 
aemlnara and the ensuing discussions will be reported more completely In the Section 3 on Aca¬ 
demic Program. The work outlined In this report was done In collaboration will Professor 
Wiener and will be presented In detail In one of the matnema - let Journals 

The theory of the beat linear prediction of a single stationary time series waa 
developed simultaneously In Russia and the United States about fifteen years ago [], 2l , 

The men chiefly responsible for this work were A. Kolmogorov and Norbert Wiener. The essen¬ 
tial foundations for the prediction of more than one time series haa only recently been made 
available and thla Is to be found In Professor Wiener's matrix factorization work [31 . Moat 
of the results reported here are essentially the matrix factorization expressed In the explicit 
language of prediction theory. 

Suppose that two slmultaneoua series or measurements are made at Integer polntn In 
time and auppose that these two series sre the measured values of two parameters, respectively, 
belonging to the ssme physical byatem. Thus, as lime progresses and a dynamics carries the 
physical system Trom state to state, the measured values of two or the parameters of the sys¬ 
tem will change. It la the Job or the mathematician Vo predict the value or either or both 
of these parameters at some specified time In the future, given a knowledge of the'r past 
history up to and Including the present. It la. or course. Impossible to do thla prediction 
precisely and under all circumstances, and thus It la also the Job of the mathcmal lclan to 
set down those criteria for a good prediction and thoae conditions under which prediction 
becomes possible with our present state of knowledge. It will turn out mat the conditions 
are also Justified In moat caaea where we desire to put prediction Into practice. 

The state of the system at any particular time, aay the present, la aacumed to be 
a random variable and. without loaa of generality, ta assumed to be a uniformly distributed 
random variable taking values on the Interval lo. 1 ] . Thla variable will be denoted by 
It la also assumed that the system la conservative bo that the dynamics which carry the 
system from one atate,a. at time zero, to another state, a', at one unit In the ruture can 
be represented as a measure-preserving point transformation, T. that possesses an Inverse 
and la such that T(o) - o'. In general, during n-unlta of time the dynamics will carry the 
point a Into the point T n (of). 

Now let fj and r 2 be two parameters depending on the system ao that at time n 

r, - r,(T r '«). 'r 2 - r 2 (T"iz). 

If the aeries of values of the two parameters measured at time zero, one unit In the past, 
two units In the past, etc., arc written: 

rjtm). fjtT* 1 *), f | (T~ 2 <x),... 
r 2 (o). r,(T‘‘o.), r 2 (T* 2 oO — 

Then a "good" linear prediction of the value fj(T”o<)i 1 - 1, 2: m>0 will Be that linear 
combination 

i\ a U«) - XT Cr.fT-VlpWd. J) 

1 J-1,2 n=5 J " 

such that I 

(1) /[f^T 0 .) - L{"> (•)]*<« 

la a minimum. The minimum value of expression (!) la the mean-equare-error for the indicated 
prediction, where the averaging Is taken over the ensemble of all possible atatea of the sys¬ 
tem. The criterion for "goodness" Identifies aa a solution to the prediction problem that 
linear combination whose mean-aquare-orror la least. 

The assumption that T la measure-preserving la equivalent to the assumption that 
both time aeries are stationary and helpsmake the following solution possible, A second 
assumption that must be made for the following solution la equivalent to asking that the time 
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URAL 'At' I'HoOl R- 'r.Al'CH 


series do not contain anj determlnlsUc components; tha< Is. any components which could De 
predicted with zero mean-square-error. 

Wltn these assumptions Professor Wiener has shown D1 that It Is possible to expand 
either function f.(»). 1 - 1. 2. In terms of an orthogoi.a. set 



9 j(T' n ») . J ■ 

- 1 . 2 . n - 0. 1, 2. 

.... 


(q,. djT-) - 

/q,(«)q jir^d-.j 

f1 n - 0. 1 - J 

0 Otherwise, 

Where q (ot) Is 

Thus: J 

linearly dependent 

0 

on r 1 (a<), r 2 (ot). r, 

(T-W). r 2 (i'l»). r,(T' 2 *), r 2 (T' 2 m).... 




(2) 

- z 

E «j< T "*>< f i 

_-n» 

• V >■ 


>1.2 

n-0 


In other words. 

, If *)0, 



(J) 

r l (T B «() - ?} 

*>(*) + a ( (m, (a.) 



■£ 

£ q,(T* , " l Jj)(r 
z n«»ir J 

»• 


-E 




J-1.2 n-0 


where Fj^'toA) 

la linearly dependent or, tne present and the past of both time series and 

r 1 (T B ’*) - r, t " > M - o 1 (m| (w) 


is orthogonal 

lo no present and 

pasi of both series. 

This, however. Is equivalent to saying 

tnat P , (m) (oc) 

la the beat linear 

predict'.on of r,(To0 In the sense defined stove. If we *se 

the fact that 


00 


(4) 

9,(«-) -Z 

£ fjCr-'kUnCi. 

J) 


J-1.2 

n.O 



for some coefficients / (1. J). a tittle computation Immediately shows tnat 


(5) 


F 


(" 


the best linear predlc- 


"’<*> - Z £ r j (T ' n *> Z I Ev ,J ' ; ”/“« (! ’ J) ) 

>i ,2 n=C5 /-.0 -'-1.2 

where o +m (l, •>) - O',. . The coefficients P n ^ m) (l, J) In 

tlon or f^T™.*) arc found to be: ^ 

(6) p n ( ”>(t..-) - Z f Z v (J ’ 

>-0 1 J-1,2 

ine coefficients I t 1. 11 are Immediately round tn terms of « (1, J) Trom (A) and the orthogo¬ 
nality relations or q,(T n «) . It remains to determine » (1. j) from some measuraole quanti¬ 
ties. 

Some direct calculations from Wiener's scries expressions for q,(«) give a (l, J) 


In series form with the terms oC the series Involving 


fc / 


and 


2 v 


V 0) %<*> 

ir 

J t.( 0 )e‘ 1 "°da. 


absolute squared value. The las'. tw<. sets of quan l'les ca i at. wn • . u. t- a d •> • >.. 
auto ar.d croaa-correlations of me oilgli.ai aeries In 'L» ! , o«* *»g «/» 

w 

ft /^ j(8).- ,n8 do - (r,T r -. r,) 


r. / J f'.(0»|*’e : "°d0 . (fjT'". f, ). 

-W 

Finally, If we assume tr.a* T la rrgodlc, ii> a- a. l at,a •.•rc-oa-cwrrt- .h ns can br computed 
Bo.vly wit n a Knowledge of the preacn* and t.;.«»• of r. origins * ’**• «.•#: 

(T) (f,T-. rj j - ^ 

> -u 

The solution for tne prediction of Hum and more time aeries generalizes 
Immedla'ely from the aoove results. 


Bayard Ranwln 

References: 

11 ] Norbert Wiener. Extrapolation, Interpolation, and Smoothing of uiatlonary Imc Series, 
The Technology Press of M.I.T. and John Wiley and Sons. Inc.. New Yor. ( 19 1 *>>) . 

[2] A. N. Kolmogorov, Bull. Acad. del. URST, Ser. Ma'n. 5 (1901) pp. 3-|A. 

[5] Norbert Wiener. Comments:-!! Rath'mallei Heivellol, 29 (195p! pp. 97-111. 


COUNriNU STRUCTURES OP FINITE K- ATI0N1. 


A Whirlwind I program for counting J.,ad!c relations on a finite set haa been written 
and tested [lj. The computations should be completed shortly. 

Besides counting simply a 1 the non-Isomorphic s' met. .res or pairwise relations be¬ 
tween objects of a set, formulas are alao being evalua-ed for the numbers of non-laomor-phlc 
structures In the following categories: 

1. Reflexive relations. IT the relation lo specified In matrix form. a.. - 0, 1 
depending on whether element 1 of the set bea-s the relationship to element J. the J 
reflexive relation has l’a on the diagonal. By writing 0'a Tor l'e, a reflexive relation 
becomes lrrcflexive, Which may be Interpreted as tne number of directed graphs (wph at moat 
one arrow Pj-ePj) on n nodes p,...p n - 

2. Symmetric relations: a ( , - a,,. 

3. Symmetric and lrreflexlve relations: a. - a,,, a . - 0. This formula gives 
tne number of non-dlrected graphs on n nodes with a- moat orto connection p^eePj. 

A. Assymetrlc or antisymmetric relations: Here a . - 1 Implies a , - 0, but a 
pair a,. . a - 0 Is admissible. For assymetrlc rolaiions 8”, - 1. for snt‘Symmetric, a.. - 0. 
Tnese formula# are given by Davis [2], 

Referentes: 

[ll Machine Methods of Computation and Numerical Analysis. Quarterly Progress Report No. 15, 
p. 10, March (1955). 

[2] R. L. Davis, Proc. Am. Math. Soc. A. A86 (1953). 


M. Douglas Mcllroy 


wiiciw 7^(5) i.aa 


pi oi»<-: -j L..a 


1L cau U dirLJixl! 


p -4, S'* 


th* o*e $eri4g css*, its 
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MULTIPLE SCATTERING OP HAVES PROM A SPATIAL ARRAY OP SPHERICAL SCATTEREHS 


General Formulation - The task of devising an electronic wave function In a periodic lattice 
la primarily that of Joining amoothly a central-force aolutlon, appropriate for the almoat 
spheric* lly-aymmetrlc atomic lnterlora. to a plane-wave aolutlon, appropriate for the almoat 
force-free Interatomic reglona. The work of Slater [ 1 ], Korrlnga [2], and others pj Indi¬ 
cates that Whenever It la allowable to simplify the problem by Baying that lnalde a sphere of 
radius a the effective potential la spherically symmetric and outside this sphere, out to 
the lattice cell boundaries, the effective potential la aero, then a technique of Joining at 
the spherlaal surfaces r . a will produce the appropriate aolutlon. Several forma of Joining 
solutions can be derived: the Important thing Is to find a form which can be easily manipulated 
for calculation. The present paper dlacuaaea a modification or the forma given earlier, which 
appears quite promising. 

To assume that the effective potential la spherically symmetric lnalde the sphere 
r • a and la tero In the Interspaces la to make the problem Identical with the acoustical 
problem or the transmission of sound waves through a regular array or spherical scstterera. 

In both cases we need not know the details or the wave solution inside the spheres In order 
to determine the solution In the Interspace; all that Is needed Is the value of the logarith¬ 
mic derivatives, at the surface r - a, of the Interior radial factors for each spherical 
harmonic. In other words, If the wave function at the Surface r - a la 

*lnt-£ (*•■•> 

wnere Y 1 are tne usual spherics, harmonics A,then all we need to know about the Interior, 
r<a, are the values of the varjoua g*a, defined In the equation 

( 2 ) 0 T 1 nt / * r > ’ * (r - a) 

m£ 

Knowing values of ihe g'a, the problem Is essentially one of finding allowed values of the 
P's for the required lattice of spheres, subject to the requirement that, as one goes from the 
zero'th cell to the p'th cell In the lattice, the wave la multiplied by the factor exp(lK.R ), 
where K Is the wave number vector characterizing the particular solution and R la the dTs^ 
placement vector from the zero'th to the p'th cell. Knowing the F's It is, oTPcourse, easy 
then to find the wave aolutlon Inside the spheres. If this la needed. Equation ( 2 ) Is also 
the appropriate form for acoustics: waves in a lattice of spherical scatterers; In the 
acoustical case the g's may be complex numbers. 

The solution In the field-free Interspace can, of course, be expressed In terms of 
plane waves 

(3) Y(r) > 21 (outside spheres) 

where K„ la the vector from the origin to the j>*th point in the reciprocal lattice, and the 
Joining equations can be set up in terms of the B's Instead of the F's. However, in actual 
practice aeries (3) converges much more slowly than series (l) so It Is usually better to 
solve first Tor the F's, then for the B's. 

The equation Tor the F's In terms of the g's la most easily obtained by use of the 
Green's function appropriate for the lattice, corresponding to a point source In each cell, 
at the point r + R , the source In the zero'th cel! having unit amplitude and tnet In the 
P'th cell having amplitude exp(lK.R p ): 

(l,) °K.k (l; k ) ’ £l exp(l S-5p ♦ **|£-So-"pl >3 /|£-Eo-" p | 

P 

- [e*P l(£p+K)-(E-£o)]/ [|* +£| 2 ->‘ Z ] 

where v !a the volume of one cell. These are, of course, conditionally convergent aeries; 
justification Tor their use will lie In the convergence of the Integrate Involving them. 

The parameter k measures the energy of the electron or the frequency of the sound 
“■•Jig ,tn ®g "* 8S91mo lhal ln -be field-free Interspace the homogeneous equation Tor Y or G 
,s v * + k t • 0. The Joining equations for r - a must, of course, give us the relationship 
between K and k, ln addition to determining the P's. 

The Integral equation for Y ln the Interspace between spheres la 

»<r> - U/*»)J^ k (r|r > )(3r/3h > ) - dS 


over the central spherical Surrace and over the cell surface. Because or the periodicity of 
aurrace <1UC, ^ ' nt ** r * 1 over the c * 11 aurface vanishes so the - o.talde r - a and or. 1 he 

»<r) - (,W)J f»<«o)(*W* 3 ) . O^IrlaJlSY »•„)] d^l 0 

where dCl Is the element of aolld angle slnfl d$d* and a la of eng- a, of direct!..'. 
g.v,-n by J . Setting r - a leads to an ln¥egPa:°equatT 8 n r.: Y or. ■ n- surra •• or the 

Sphere, which Ta the equation we desire. Using Equations ( 1 ) and ( 2 ). we obtain 

• •„£, S*v*w 

+ andn °K.f 

which la the fundamental homogeneoua equation determining the F's and k ln terms or K. 

Since Equation ( 6 ) la a homogeneous integral equation, It nos non-zero solutions 

for only certain discrete valuta of k for a given K. These sre the 0 owed .evcla we seek. 
»e can find this relationship by a variations] mctliod, as Kohn and Restorer [ 3 ] suggest, or 
else we can solve the Implicit set of simultaneous Equations ( 6 ) directly Tor k, for any 
chosen K, since the P^'a converge quite rapidly. 

To put Equation ( 6 ) ln useable form, we apply the uaua [3, k] expansion of the 
Oreen'a function (for a<r<R ) 1 

P 


Ogx^lfio) -.(1KA») 5Z c'lLSp/e'AlVVll’dn „ 

— p B 

- u£ (81+O^f Y?($. m){-Yj(5 0 .m 0 ) 


z 

ml 

* E ! 

m'/' 


J|(ka)h^(kr) 


l-l 




23 £ 'n'"'< 0 p-V'p)e'!i-L p h n (kB p )} 

where 0 , 0 are the spherical aisles Tor the vector R and where the I's are the Gaunt factors, 
deflned p by The equations ~ p 


(7) 


( 8 ) 


Vj(ff.m) #($,») - 23 I n (mf |m'Z')Y”'-(^.m) 
n 

We also need the g factors for the Spherical Bessel functions 

SjUa) . -[a/Jj(ta)] [dj f (ka)/ds] - tanp^lkajj 


These would be the.values of the g'a for Equation (2) If the effective potential Inside r 
were zero. Tables of cW^ are available [ 5 ]. 


(9) 


Substituting all this ln Equation ( 6 ) we obtain 
JF m / - wa J i (ka) [Ugj-gjJh^kalF^ 


m'/‘ 


£ .^'-(^elS^ n n (kR p )]F„ 1/ , 


PWO 

where, for each choice of K, k la adjusted to make J • 1. This la a compact and rapidly con¬ 
vergent form for uac ln computations. The quantities ln square Brackets do not depend on the 
spherical lnterlora and me be computed, or.ee for all. for each type of lattice. The spherical 
Interior only cornea In wltn the terms !g,-g5). wr.loh arc, as wc one 1 see later, aloaeiy re¬ 
lated to phase angles for scattering fro# a single sphere: These go rapidly to zero for In¬ 
creasing l ar.d ln many cases only I - 0,1,2 s Trice for throe or Tour place computation of k 
and of the F'a. The flrat-order effecta of higher I'n may be added, if desired. 
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It Kill oe shown later that for regular lattices the matrix or the coefficients of 
the P's Is He ml lean as long as the g'a are rea , so tne compiled values of J are real when 
end K are rea . We will also show that When x +K the quantity In sq are brac.eta becomes 
Inversely proper-tonal to K -x . so that whan the scatterer Inside r - a la wea< (or when 
a-» 0 ) and a the g's approach the g 0, s. then approaches K, as It should. 

Furthermore, the Equations ( 9 ) are exact, as long as tne initial assumptions 
(V - 0 Tor r>a) are valid. Although the derivation given stove assumes a reg..ar lattice, 
the same equation can 00 obtained without the assumption of regularity, as ong as the mean 
density of spherical center* (l.e.. the mean density of the points represented By tne vectors 
H ) approaches a constant for large R. The summal Ion over p la then the s imitation over all 
The other scattering centers, ’he p'th one being at R . Effects of lattice distortion may 
thua be computed from the same equation. 

Prom this point of view Equation ( 9 ) la an extension of tne method of Poldy [h] 
rot' calculating the multiple scattering or sound from s random distribution of scatterers. 

The quantities In the square QracKets represent tne effect, on the surface of the central 
sphere, of the scattering from all the other spheres. 

T.. Q -o ttlce Sc attering Function - The Basic lattice scattering functions [ 7 ] M* are defined 
as follows: 

(10) t 1 * 1 ITe'S-Sp V> p .0 p )h n (KR p ) - -1 S 0 J 0n * 

It turns out that as < aoproaches K, M goes to Inflnlt; as («»/ku) (-K K) n Y"j(u , it )/(K 2 -k 2 ), 
where 1 / Is tne cell volume and u and U are the apherlcai angles for 

vector K. In many parts of the calculations It la advisable to subtract It out. The func¬ 
tions My as defined, at- real quantities for reg.ar a '.Ices; they have a sequence of poles 
when x*- •|K+K i/ |2 and except near It-: poles tl.e. are fairly well-behaved f-notions of x, 
having magnitudes of the order or unity. 


of the order or unity. 


The method of calculation or M Is, or course, a generalization of the one used by 
Ewsld [B], else by others [2, 5], We a.bslt.te the source distribution 

e<?> ■ f" i;(2V n+5 /.i , ')|r-R p | n Y-(J p ,a p ). 

p^Q 

• exp [tK.Rp + (ic/Zp) 2 -^IP-SpI* 

(where .9 , » are the apnerlcal angles for the vector p-H ) Instead of the set of point 
sources,RpteS). which prod.ee the aeries on the U-rt-nan3 p slde of Equation (10). Adding the 
term P - 0 to p(9) a. lews the waves produced a- he origin by this distribution to be ex¬ 
panded as a Fourier aeries, giving for part of M 

<>» C - *-»>■ ft fEfUJWf ^0^-^ flev V) 

' 3^*7 *>•' Vo) ( 

where, as before. K-, Is the vector distance to theb"th point In the reciprocal lattice, u„, 
“p are the Spherical angles for the vector K + K,, and v la the volume of one cell. Wc have 
subtracted off the pole at k‘ - K.‘ for ease In computation; It will be added separately. 

The wave from the p - 0 distributed source must then be subtracted off. Thla la 
zero except for m - n - 0; the aecond part of M la 

(12> t * <U„[‘ * 

- z -»>] 

- KA* I‘A*) f r«V*!^ 2 , ' l /al(*a-l)] 

8a O 

b *iJ e u n l !j e w,,ves ' ac r ■ °- djc to th « «ctual P(.o) and to the aub- 
stltuted has to be added. The p'th term of this correction series Is 


OKADUA h aCHOOl HMKAflCM 

( 2 -V' i - n )>X(» p .fl p )«pf'K.P,. • 1 2 ^,] 

S n-^ 5 /fV“ R P »%«“ r > ' "n<'\>-n« V*V- 
Rp - 

- ta K( r x p .^r p ) 

- (;x)]ox 

x p 

nSSctra; T# -“ ,r, ° wr 10 * - v - - 

S n =e'*p 5(2pP. p )- <r in{h n (RR p >F> p) ) 
where F^a) - <W, . * 

Prom thla we finally obtain for the third part of K 

(15) " 2 n~ y 2 Yl ( f ft p ) n * 2 w,p<lK.K p - /r 2 ). 

P P " P 

If 0 is as large aa 2 divided by the distance to nearest neighbor, thla sum converges very 
rapidly, and usually only • he nearea- nelgnli.t terms nesd tt.- alned if four-flgun 
accuracy la required. 

Aa an example, for a simple cublo lattice of spacing 2b. if we set/?- (tt/B^b) 
and use the dlmenalonloas parameters J 

z - ibK/r); Z - (bK,tr); R^ . R p 2b) - n x i + n y J + n^W; 

Zj, - (bKp'wl - Xi ♦ lex(h x . n y . n ix . A . & Integers) we have 

(i*) (z.z)^ i z**BL - i im4l +2*b 

° A ? z V \l + lj - * 2 

+ (2^ m^z)exp(2z 2 - gm S )[cos(2»Z > ) * eoatevZy) 

♦ coa(2rZ )] (1 - J±S«!) + s* V 

tr 2 a ! ( 2 s- 1 } 

8-0 

The first series requires about a hundred terms to obtain four-figure accuracy, but the task 
la not difficult to program for machine computation Curves Tor M)j are given In Figure 1. 

The polos of the function come Trout the terms In the rirat series. If the terms Tor the 
nearest two poles be subtracted from M, the remainder can easily be Interpolated to obtain 
Intermediate values. Table I shows s few representative values. We note alao that M° for 
* (£4-na x ) equals Tor Z. ° 

General Properties of the Solution - It Is now possible to simplify further Equation (9). 
we note that only for T' - l, m' - m la there a tenp In In the aum over n. Furthermore, 

I (m/|m/) - (f+m) !/( 2 /+l) (X-m)!, so that the -iJ S ter* can be combined with the first 
term in the curly brackets In (9). Since h, - J? + Tn, the J, term is cancelled out. Since, 
alao, from the formula for the Wronaklan ToT »nd n^T and rfom the definition of g£. wc have 

(15) -ka(g / -g|) J|(waln|(ica) - 1 - xa^(kalnjfka) [g^- tan ^(ka)] 

where fo(*) - tan* 1 [-znj(z)/n^(z)] 

Tablea of have been published [-]. 
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Table I Director! I. K . 5 j( Ki Simple Cubic 



Values of 

*"o° -(£){ 

1 . _ + 





Z 2 - z Z 

(1-Z)* - z 2 



z 2 

Z - 0 

0.1 

0.2 

0.3 

0.4 

0.5 

-0.4 

-.7712 

-.7758 

-.7795 

-.7823 

-.7839 

-.7845 

-0.3 

-.7110 

-.7164 

-.7206 

-.7238 

-.7258 

-•7265 

- 0.2 

-.0467 

-.6530 

-.6580 

-.6617 

-.6641 

-.6649 

- 0.1 

-.5775 

-.5850 

-.5909 

-.5954 

-.5982 

-.5992 

0 

-.5023 

-.5113 

-.5185 

-•5239 

-•5273 

-.5285 

+ 0.1 

-.4196 

-.4306 

-.“395 

-.4462 

-.4505 

-.*519 

0.2 

-.3270 

-.34 09 

-.3520 

-.3607 

-.3663 

-.3b82 

0.3 

-.2214 

-.2393 

-.2537 

-.2651 

-.2725 

-.2752 

0.4 

-.0971 

-.1211 

-.1405 

-.1561 

-.1665 

-. 1702 

0.5 

+.0457 

+.0124 

-.0148 

-.0372 

-.0526 

-.0583 


Values or 

+ pJ z 

1 - Z 1 




1 W/Z ? - z i 

0-Z) 2 - z 2 J 



-0.4 

-.03619 

-.02944 

-.02249 

-.01524 

-.00771 

0 

-0.3 

3897 

3168 

2422 

1644 

833 

0 

-0.2 

4222 

3429 

2625 

1785 

906 

0 

-0.1 

4606 

3737 

2867 

1956 

995 

0 

0 

-.05066 

-.04102 

-.03151 

-.02156 

-.01099 

0 

+0.1 

5629 

4546 

3503 

2407 

1232 

0 

0.2 

6333 

5097 

3942 

2726 

1401 

0 

0.3 

7237 

5797 

4508 

3142 

1625 

0 

0.4 

8444 

6714 

5263 

3708 

1932 

0 

0.5 

-.10132 

-.07968 

-.06319 

-.04517 

-.02376 

0 


Values of 


_zf_ 

, (i-z ) 2 

1 



[a* - 

(1-Z) Z - z' 

'1 


-0.4 

-.03619 

-.03416 

-.03226 

-.03065 

-.02955 

-.02953 

- 0.3 

3897 

3711 

3521 

3349 

3227 

3230 

-0.2 

4222 

4061 

3869 

3682 

3544 

3552 

-0.1 

4606 

4479 

4286 

4077 

3916 

3930 

0 

-.05066 

-.04991 

-.04810 

-.04553 

-.04359 

-.04378 

+0.1 

5629 

5628 

5426 

5136 

4892 

4915 

0.2 

6333 

6443 

622e 

5865 

5547 

5570 

0.3 

7237 

7518 

7278 

6798 

6365 

6381 

0.4 

8444 

8996 

8704 

8031 

7415 

7406 

0.5 

-.10132 

-.11070 

-.10712 

-.09747 

-.08846 

-.O 8787 
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Z - (bK/ir); z - (bn/ir). 
Lattice spacing 2b. 
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( 10 ) 


Finally Equation (9) Becomes 

[(J-l) - ^aJjtKalnjfKaltgi- •"/V] K m/ 

+ RaJ|(«a) * * 6 t >" e i >' •*! • lua) 

m'i* 


[ (2<+1 > {zmH S l„(m<|m'*')M^ - 0 


where we adjust trie vaxue of < ao that J-l vanishes. The Hermltean character of the matrix 
components Is apparent. 

The factors In g are closely related to tne pnaae angles for scattering from single 
spheres, for 

(17) 


tanJJ 


J l (^a)(e e -gf) g f J/(«a) + ka.l^Us) 

nJOaTrgj-'.an jj " zjt'.faa) ♦ "•'•nJPsT 

In fact, when all then's are small except for £ , we can approximately reduce Equations (16) 
to the one for l. m m 0 and, setting J-l -0. ofcraln 

k col £ t*kM° 
o o 

a an Implicit equation to determine for a given va.je of K. To thle approximation, 
re, the electronic enprgy k2 (or the mean Index of refraction K k for ao und waves) la 


(18) 

which 1 

therefore, the electronic enp 

determined by a term k not S .. depending on the Interiors of the lr,dividual spheres, and on 
a term M®, depending on the lattice geometry. 

« Por energies Just above tne patent!*, of trie Interspaces [t 2 sms 11 and positive) 
k cot J varies o.owly with K or k and H? may be written as (li/v) (K 2 -.: 2 )* 1 + N“ wncre NO 
la an,aIt. The form of the equation for“lteratIve solution la 

k 2 - K 2 - (u r v)(k cot cf - N®)' 1 

Near a resonance level of the atom .< cot S goes o infinity: In this caae the rorm ( 18 ) la 
rror>- appropriate for- Iterative acl '.on. near the edge or a Br'.llouln zone M° Becomes very 
large, no that this la the dominant term. Having curves for k col <f and kM° as functions of 
for s given K, graphical solutions car. ne obtained fairly easily, 8s long 8s only S la 
largo enough to be worth Including. ° 

When J, for t >0 Is small but no' negllglnlc, we can solve Equation ( 18 ) first to 
obtain a firs- approximation for , substitute tnls in the following equation Tor the 
coefficients P, 

,l9) P mi) “ l/(2 **"{f;-^7 ^T k^TVanffw l [ |cot Jo-"o ) l- l >*- 


* «rh 


where, aa mentioned before. , . v are t. e spherical a, g - s for the vector K. A next 
approximation for car. t .-r De Obtained b, subttltutlng Back lr. M.c eq .a'.Ton for F 


When ,! 'r f are not so small :*.»• 'heir squares can he neglected, then the determl 

nerr. ol the coolr 1 r 1 cuts ol i-.q.atloa ( n) must be Bet equal to zero and solved In Its entirety 
31r.ee, In moat can. a of Interest, gi can be neglected for f?2. It Is seldom neoeaeary to 
solve a determinant of .arger than third order. In practice, the approximations of Equatlona 
1 ■ 1 a "'> < •"» d'-l" good r • •«. elements of t. fr-sr r w of he periodic table fol. end 

a two-by-two determinant la Sufficient Tor sodium to potassium. * 1 

It la sometimes convenient to expand the solutions In the field-Tree lnter-apsce In 

P'ven K E ,, ;: 1 :: V S! Tr i8 ,'r T " "'•» •«» - M*t been compute Carl 

A u K ,'„ l t .? '? '■■■ p "«ve vxpsnslon for 0 given lr, Eqrstlon (A). The 

aecond term In he Integrated ol (5) ira. eaal: be Integrated, -sing ( 2 ) and a rein* ion 1 A 1 
neiweon he Vs and -a. The firs, ten. oa. .1.0 be .l-pllfi^^li; ttewtatlS 11 

aos, ‘V aa ^o* »lnu p air, cosf^-v^,) e’fK+K^).^ 

- *lrl(-l) X Jj(|K+K,|.) Y"<^ >v ) 


inablnlng the two, we see that. In the Interspaces 

(TO) » 


-.KIM,.,]. 

• J/( | 5 +!i| •)V?(u y .v t ,)K r } 


This equation can, of course, also be used as a basis [ 3 ] for a set of eq .stlons 
determining the P's or the coefficients by of the plane waves or Eq.ua'ion ( 3 ). Bu- ... con¬ 
vergence of this series Is much leas rapid than the series In the P's. 0 It '.a much easier 
10 use Equation (lo) or ( 9 ) to determine k arid the F's and then compile series 120) If needed, 
rather than to use Equation (20) to determine k and the P'a. 

Tables of the It's are now being computed by the Project. 

Philip M. Horae 
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[a] P. P. Ewald, Phys. Rev. 62, 107 (19*3)- 
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ENERGY BANDS IN GRAPHITE 

In an earlier report [l) an unsatisfactory tight-binding calculation of a two- 
dimensional graphite model ijas described wherein the overlap matrices were not positive defi¬ 
nite for certain values of k, the reciprocal lattice vector. A tentative conclusion was 
drawn that the cause of the difficulty was a result of only Including overlap Integrals up 
to third neareat-nelghbora. Since then, considerable time baa been spent writing and testing 
a more elaborate matrix generation computer program, which allows the Inclusion of up to ninth 
nearest-neighbor Integra.a In both tne Hamiltonian and overlap matrices. On the basis of the 
new matrix generation program, the hypothesis or Ina-fflclent overlap Integrals was found to 
be the correct Interpretation of the previous unsatisfactory calculation. Further Investiga¬ 
tion la being made to determine the number of overlap Integra a which are significant. 

An additional tea' hea also been made to Investigate the Importance of the hereto¬ 
fore neglected three-center Integra,* or the Hamiltonian mairlcea Two seta or results, con¬ 
sisting or the energy eigenvalues st severs! representative k vaiuoa. were obtained: The 
firs' set was computed with all three-center Integrals neglected; the second set was computed 
with estimates made or the three-center Integrals. These 'wo sets of results differed markedly, 
e .d therefore It appears that at least the most Important three-center Integrals must be 
accurately Included in this type of calculation. 

At the present time these required three-cen'er integrals are being computed using 
the multi-center Integral techniques desc-lbed briefly In tne previous progress report; namely, 
that or expanding all wave functions and potentials about a single center and performing the 
resultant one-center Integrals numerically. Pina testa of the computer programs necessary 
for this procedure are currently being made. Preliminary results indicate excellent agreement 
with the previous two-center Integra, values (obtained by the different technique of evaluating 
a.x.jtlo functions). Thus, the extra .sbor or obtaining me uirea-canic. Integrals has some 
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oompcnaatIona In thgt It rorms an Independent crosa-cheon of a large part or the numerical 

work. 


Fernando J. Corbato 


Rorerenc®: 

Ill P. J. Corba-o, Quarterly Progress Report, SOUd-State and Molecular Theory Oroup, 
October 15. 195*. P- 5. 


VARIATIONAL DETERMINATION OP ATOMIC WAVE PUNCTIONS 

The Whirlwind program for carrying out the variational calculatlona for the florae- 
Voung three parameter ware functlona hao been nearly completed. A thorough recalculation of 
tne original energy tablea aa well as the calculation of aeveral new atatea such aa 
(Ia22a22p5)2p will be started shortly. 

The gradient method mentioned In Progress Report No. 15 works very well as long as 
the parameters are of approximately equal strength. This condition la not satisfied for atatea 
with many 2p electrons and we have a situation like the one Indicated schematically In Figure 1. 



Figure 1 

Energy Contour Diagram 


a. b, and c represent the successive points determined by the machine In trying to reach the 
extremum. The narrow hill makes the path to the top extremely Indirect and lengthy. 

To remedy this, the program has been modified In a manner illustrated In Figure 2. 
An Initial gueaa (A) of the beat parameter values Is made and the gradient la calculated. We 
then proceed along the gradient line until we reach an extremum (A 1 ). At (A') we again cal¬ 
culate the gradient and move In Its direction to (B). Starting at (B), we the n move along a 
line parallel to AA 1 until extremum (B*) la found. Finally, we move along A'B' until we 
reach (E') which la taken as the approximate minimum point. This process may be Iterated if 
necessary. The varloua minima are found by linear Interpolation of derlvatlvea obtained by- 
taking first differences. 

When applied to states with five or six 2p electrons such as (la 2 2a 2 2p®) 1 a, the 
above procedure gives satisfactory results. 



Figure 2 
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An effort la now being made to improve the wave runctlona by adding extra para¬ 
meters. The moat attention la being given to the la and 2p functions. The forma now being 
investigated are 

1. . Mj (.-/•*•♦ *-/*►) 

2p - N 3 (B'rV'‘ ; ’°»-)pJ(oosf) 

The values or B and f have been calculated Tor the helium lsoeleotronlc sequence and It la 
Hoped that theae will hold for any configuration containing two is orbitals. 


Arnold Tubla 


SCATTERING OP SLOW ELECTRONS PROM ATOMIC OXYOEN 


The physical aopccta of this problem have been discussed in tne various theala 
reports. Roughly they can be summed up aa an extension of the Hartree-FocR method to a con¬ 
tinuous state problem. Mathematically the problem la to solve a set of aeoond-order lntegro- 
dlfferentlal equatlona. The equations are of the form 

i-jj + r(x )u(x) + g(x, / p(x')u(x')dx') - 0 

dx o 

p(x') la a known function. The constants C^ are Integrals from 0 to ao over tho unknown 
function: K 

J x 1 ' p(x' )u(x' )dx' 

o 

(kl 

To solve theae equations, the procedure Is to take a guess at C ‘ and to solve the equatlona 
using these guessed values for the c'*). This gives the function u fro® which one can evalu¬ 
ate tne Integral Cl*) and use this as a next guess; l.e.. one uses an iterative procedure 
which presumably will converge. 

The equations solved thus far were s-wave equations for the total S » 1/2 case 
(spin of Incoming electron opposite to the S - 1 spin of atom). In this case the angular 
Integrations make the Integral part of the equation vanish and the equation la Juat an ordi¬ 
nary differential equation. 

u {<a o R) Z .- SOjU)}- 0 

G(x) la a known, tabulated function which depends on the numerical Harr.ree-Pock wave function 
for the oxygen atom; (a k)2 is proportional to the energy. Actually we are Interested In the 
phase shifts and these- have been computed by- Whirlwind to be. for various values of the energy 


(in radians) 

(»o"> 

Energy (in e.v.) 

-5.11 

.1 

.136 

-3.09 

.2 

.5*5 

-3-07 

• 3 

1.22 

-3-07 

• A 

2.18 

-3.07 

.5 

3.*o 

-3-08 

.6 

*.90 


Aaron Temkln 
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NEUTRON-DKUTRON SCATTERING EIGENVALUES FOH A SPHEROIDAL WELL 


I-. ...a been found ir.t In- r.-V ici::»r:.< .-ngma are much .or. aermltlvetoth- 
aeeuraey of the ass ,med deuteron wave function tMi wee previously believed to be the case. 
P-um Kq.a-lor.a I . 0 ) ar.il ( . ») of tile previous report [.]. we can for™ the Integral express.on 


M, -f*s 9jTt«k-<V 


II 








S 4} 


■ 1 . 011 - 01 . to,, S . a, ' I VO. ye . do.-- f .net Ic aid ar defined In ( [.], 15. lfl). On 

appl. ng Ore.-n'S theorem r., a--: i.d de:-'va ,, ves e-.u.’.n 1 ned 'tv T, we can write (1) In the 
form' 



■ n.-re L. la the surface Integra defined In It.] 27o) ar.d la another surface Integral of 
the Bar" type. Equa Ion ( 2 ) la ae laflrd Idem lea .. if uj_, w . w,, we use the exact 

deuteron function; »lt: a approx 1 ma• ■ d-.-uter n function. in* left and rTght-hand aides of 
the eq.alor will differ. Th'.a difference leads ua to trie font, iat tori of a variational 
principle slightly different 1 'rom trie one dev-.op-d In ( [.) 2 -u). which ■•■ ana. I denote aa 
A. The alternate form, which we denote aa h. la obtained by adding ( 2 ) to eq.atlon ([ij26) 
and then proceeding xact ly aa In the origins derivation. Slice tne terms added do not tn- 
v v- the trial I merlon the a' at lonar 1 of the resulting «q .at Ion with respect to varia¬ 

tions In y la unaffected. Varlat tonal principles A and B are, of course. Identical In the 
limit of exact, deuteron function. 


The deuteron wave function which has been raed this fan In the calculations la of 

the form 

(3) F(r) - i* rZ + ce* r " 

when constants a,a, c are determined by ar. Independent variational calcu.atlon to mlnl- 

mlz" ... of the' deuteron ground ata'e. (A Oa .aalan form for the deuteron function la 

required In order to perform the Integral Iona.) Tr.e form (3) gives a finding energy which 
differs from the experimental va.-.e i- ar.out .2 MEV, which la quite small In comparison with 
the kinetic mid putentla energies Involved. It was therefore be.leved tlrat (3) would be 
• n adequate function 1 o .ae lr the rta lev. la'Ions 

However . in j.aUon ( 2 ). wt icr. la not a'a’lonary. the uec of ( 3 ) as a deuteron 
function leada to o difference or about 20 per cent between the left and right-hand aides. 
The reaaon for this large difference la mainly that I— turns out to Involve the correlation 
between the deuteron function at smal aid large values of the argument. At large distances 
the wave runotlon (3) la the poorea'. a fa ' which doea not greatly arrecl the binding 
energy . 


Aa a result of this 20 per cent difference In Equation ( 2 ). we find that the 
variational principles A and H give results which als< differ Substantially. Various plausi¬ 
bility arguments .ead .a uo believe that A la more accurate, but are not completely con¬ 
vincing. In an erfort to resolve tne uncertainty. It has been decided to redo the calculations 
.slug a mote accurate deuteron function. Christian and Oammel 12 ) give a deuteron functltn 
of the form 

(*> P(r) - e-* r2 ♦ e l .^ p * ♦ 

which approximates the exact aolutlon within three per cent and leads to a difference of only 
five per cent between the two aldea of Equation ( 2 ) (4) should therefore be completely ade¬ 

quate for our purposes. 

The use of (4) aa a deuteron wave function Introduces no Intrinsic complications In 
the computation. It doea. however, Increase the number of terms to be evaluated and will re¬ 
quire additional programming and machine time. However, the relative eaae with which the 
stationary point was found using wave function ( 3 ) Indicates that the additional machine tlare 
should not be excessive. The reprogramming la now In progress. 


References: Leo Sartorl 

[1] Machine Methods or Computation. Quarterly Progress Report No. 15 . p. 23 . March (1955). 

[2] R. S. Christian and J. L. Oammel. Phya. Rev Jl. 100 ( 195 , 3 ). 
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In order to compute the matrix elements deaer'.oed lr. the 1 -,-v > .s q ,ar - , , -pet • . 

It la necessary to compute the logarithmic derlva' Ives of • r,- » "rad'.a " a;.. la fur.c 

tlona Je af (h,a) he^(lg.a). Aa before, we have 

h 2 ♦ g 2 - H 2 : H 2 - d 2 V o 

where d la the rocal length of the (prolate) spheroid and v la tne w . depit, while a Is 
the value or the radial parameter at the bo .r.dary or tr.e ■»■(.. Th. - , runction is I In' e 
at a - + 1 and behaves asymptotically aa a spherical beaael fur.ct Icn. Tne .(!(,.a) blows 

up at a - 1 and has the asymptotic behaviour or a spheric*. Ha.-ixe. fund it. f ' r,e first kind 

and of argument Iga. 

Now It can be shown tna- me Je « function may be expanded lr. a aeries or spherical 

Beaael functions aa [1] t 

(1) ' X* ' 51 *x.*pf h M 

p • G for (/-m) even *~ c 
p • 1 for (/-m) odd 

Tire a 'a are computed on a demand bat la by a Subroutine which was developed b.v Mr F. J. 
CorbaCo and Dr. J. D. C. Little, Computation or the spherics Beaael functions was carried 
out by the method described In part II of this report. 


By examining the differential equation Tor the "radial” function [2]. one sees that 
tne prolate spheroidal r.nctloii he .(ig.a) la Identlca. '0 the ol la'- spheroidal function 
he .(g. la). Although this funetlorrmay e.so be expanded in a tier lea Birr I lot ,u iliar In 
Equa'Ion jl) (with the J 'a replaced ly h 'a or the Tlrat kind). It turns u.t that the series 
converge very slowly or Rot a' all. Consequently, It waa neceaaary to find a different method 
for computing this function. 


A suggestion aa to a possible method may be obtained from the asymptotic behavior 
of the he . function. Suppose we choose wo arbitrary numbers OLc-Ami ant * a va ^ ue such 
that zx>B1 Now let » , be the value and P. th. derivative or absolution of the radial 
tieroldal eq ration at"* . Then wo can wf'Te 

“-/(^1 5 * iRrfuyl J < ' ■ r <4 *• •*/J 1 J <' *•) 


apher 


?>., <*.> * ( */■**.) Li,t'**,'-!,-/.) •r'lA-X 


Now ne 
z. but 1 


.(g,iz) not only behavea as a decreasing exponential for sufficiently large valuea of 
11 must be a monatonlc decreasing function throughout Ha range. On the other hand, 
jc ,(g. Iz) must increase monotonlcaily throughout lta range, going asymptotically aa eV^fx . 
TnW. by choosing z sufficiently .Urge and Integra- '.ng the dlfferen- '.a equation for the 
radial equation lnwtfrd from z to a. ta.-lng piL/(* 0 ). P_g(* 0 ) aa lnlt,a conditions, we can get 
to a high degree of accuracy That 

Pw.,U> t i,t 3 ; xKx‘1.'"Jx.A 

Moreover, since It la only the logarithmic derivative of he^j that la of lntereat. the value 
of tne constant a la of no concern. 


c» 

w 

A* 


This scheme la now In use for computing the desired runotlon. The constants „ 

,(z ) are ao chosen that their ratio />_f(*„)/«Cx(*J * -8 which guarantees that a /a, 

*- - - - .-/IM—r°.—M—0.- s , 1 th an increment 

waa found that 


This acneme is now in use .or rump.Lina u..o ura.ice 
>,(*-). A mf (* 0 ) »re ao chosen that their ratio- -6 which 
-Hi Be arsBll to begin with By using a accond o«or lntMiaTlon process 
,z determined ao that gAz • .025 and choosing z 0 ao tha' g(* Q -«) 2 S, It w 

could be obtained with an error or leaa than about three per cent In a time or about riiteen 
seconds. 

II. Computation or Spherical Beaae; *unttlona‘ [2] 

The spherical Beaaol. Hanxel. and Neumann funotlona of order 11 all obey the equation 

(» U i, - Li- n - 


1 Equations (!)-(•») may be found In [3. pp. 1573-b] 
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or. equivalently, the difference equation 

in particular ee aha: take the two Independent aolutlona of (l) to be the 

function. J„(x) and the spherical Neumann function n (*)- These functions oay he conveniently 
defined In Senna or Inltla. conditions on he dllTerPnoe Equal.on (ZJ. .nus. 

( 3 ) j.j r*> - e <p‘ 

j o (*1 = %~ 

In addition to Equations (l) and (2). we snail also mane use of the relation 

(*) /**"*! U ,, '' x * ix> 

If we act y n (x) - j* n (*1. Equation (l) becomes 


{ fi. * ii- l’ 

Thua, we aee that the n-x plane may be divided Into two regions according to 
(I) x 2 * n(r*l); (II) x 2 >n(rv+l) 

and the beliavloar or the aolutlona will be *aaentlal> different In each region. In fact, it 
la not dirricult to ahow that In region (I) ij (x) la a monotonlc Increasing (poaltlve) func¬ 
tion of x and 1/x n (x) la a monotonlc decreasing (in magnitude) negative function of x. On 
the other hand, botfl functions are oacillatory in region (II). 

In moat of what follows we snail be concerned with the behaviour of the two func¬ 
tions in region (I). The discussion will Involve two numbers N.t which are defined by: 

(dj) For fixed x. N la 'he least Integer Such tnai N(N*l)ix Z . 

(dj) For fixed n. t la any positive number such that t Z < n(in-l). 

Ualng the above definitions together with Equation (*) and me known behaviour of 
1/x J_(x) In region VI). we see that 


1 /x J n (x) in region 


we see that 




which leads to the relation 

(1.1 


A trivial modification of the above argument leads to the conclusion that 

( 3 «) "«<>> < x"-" (ill ?; r ' 

( x) 1 W -• 5 / • 

Tn order to raa.<e use of theae lnequalltlea, we define for n>?N 
*/>*.,,<*>1 
* • *[ j„ (*> * Pr. <*> J 

where J ,(x) and J (x) are any two arbitrarily chosen numbers. Repeated application of 
Equation *2) then gives 




Now^make the apocla choloe J^jlx) - 0. Then, by making use or the Inequalities (5). we 

(6) / < X*"— 

Table I gives a few values or n va. N with n ao ehoaen that J<- 10 ' 10 . 

Table I 
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The inequality ( 6 ) suggests a method Tor computing sets or functions j.(x) (x 
fixed, n, i km 1 which appears to be Ideally suited Tor a high-speed digits: computer. We 
first put the lirger or the two integers (n_,N) In place of H In (b) and choose r. to mako S 
leas man some specified amount. Next take‘j (x) to be any convenient number and J ,(x) - 0 
The recursion relation ( 2 ) la then used to geHerate_the set of numbers T (x) agtfitt which 
la left stored In the machine and also the numbers J (x), J (x). Uae of either Sr these 
last two quantities together with the appropriate on* of the‘Equations ( 3 ) allows ua to com¬ 
pute the normalisation constant. 1/j , and, conaequentlv. the J k (x). 

In tne event tnat x (and. In consequence, N) la very large. It will be necessary 
to generate many values which lie In region II In order to obtain J . We have seen tnat use 
of Equation ( 2 ) In region I leada to "damping" of the round-off error Introduced at each atep. 
The same effect does not occur In region II. Instead, we should expect the successive con¬ 
tributions of round-off error to add In random rssnion with a resulting loss of acouraoy In 
the computation of the normallxatlon constant ol. My experience has been that this effect la 
negligible for x as high sa 25 (computations carried out on Whirlwind ualng 2 A, 6 arithmetic). 

In carrying out the computations. _a alight modification or the above method was 
used. It was found that cases occur where J (*) cannot be ehoaen sufficiently small ao tnat 
J„(x) does not exceed the storage capacity oT the machine. Thla difficulty was avoided (and 
tne computing time decreased) by ualng Equation (2) In the form 


/*<•! - *C 




and computing with the ratios throughout region I. 


References: 


.’ack L. lire t aky 


[l] F. J. Corbato and J. D. C. Little, Machine Methods or Computation and Numerical Analysis, 
Quarterly Progress Heport No. 11. p. 37. March (195'). 

[21 c.f. J. Melxner and F. W Sonafke. Mathlenache Funktlonen and Sphaerold Functlonen. 
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[3] P. M. Morse and H. Feahbach, Metnoda of Theoretical Physics, Mc-Oraw-HlU ( 195 *). 

A SELP-CONSISTEN DETERMINATION OF THF. NUCLEAR RADIUS 


A Hartree-Poch typo of self-cons latent field problem has been set up and solved for 
a standard heavy nueloua of 23 - 2 N - A - 18 *. 

The purpose of the calculation la to determine the relation between nuclear matter 
distribution and the resulting collective nuclear potential. Aa a rirat trial, lnternuclear 
forces were considered to be Wlgner plus Majorana, the radial dependence of the force being 
Gaussian: 


*3®"’ 5 rZ [l+xP( rr’ j] 


P(rr') la the Majorana exchange operator, and x la the ratio of Majorana to Wlgner force. The 
following la an outline of the method of solution: 

(a) A square well or radius 8.2 fermla (• 1.* x 18 * 1 ') and depth 32 Mev la filled 
with 18 * non-interacting particles: *6 neutrons w"-n spin up, *6 with spin down and “6 peotona 
with spin up, *6 with spin down. 

(b) The wave functions Tor thla well are used to calculate tne density distribution 
or the particles. 

(c) The collective potential generated bj those la calculated using the assumed 
Internucleon potential. 


from these 


(d) A new set of wave functions Is derived from the new effective potential, and 
wave functions a r.ew density distribution la obtained. 

(e) Steps (c) and (d) are repeated until the potential la the same for two Iterates. 
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found that tho process converged after a »;"8le neraUon when * aaa allghUy 
The significant result of ' no oaleulatlon la that for tho 


It was 

greater than A. The algnlfloant r-a-ii. w. v..» __ ..... „r , h ,, a.inir 

solution* the meat square radlua or the potential la A 1 > ferrnl*. and that of tha matt r 
dlatrlbjtlon J8.A ferotla. ao that the roo mean square radlua Of "ifSulaflSM ollh 

i,o! neni laigor than tho matter distribution Which produced It. Further calculations with 
YuKewa forcea and other oxohango admlxtureo ore now In progress. 


Manuel Rotenberg 


w.. -m rn!" r,nUh the «haly*ls of this particular building In the next two eee.e 

We will then use the eame program to analysis a similar building. 


Charloa w. Johnson 


2.3 Final Reports 


DYNAMIC BEHAVIOR OF SHEAR WALL TESTING MACHINE 


MACHINE SOLUTION OF THE DIFFUSION EQUATION 


The shear w» i testing machine la composed of two trusses and connecting bracing. 

Tho ooncrete well to bo tested la placed on the frame between the two tmeses and the dynamic 
load la applied by a hydraulic system. 

For the elaatlo-plastlo atreea analysis the testing machine and the concrete teat 
wa.l arc represented by a twenty maaa system connected by elastic-plastic springs. The masaeB 
representing the testing machine are located at the truss Joints and are composed or 1/2 or 
the maaa of eaah member coming into the Joint. Since the truss Joints will all be assumed 
oln connected (no moment resistance), the resistance or each member will be composed of Just 
the axial force. 

The behavior of each mass point can be represented by an equation of the form: 

"!*!<*•> ‘ P XI< l n> - £"xi<‘n> 


where: 


Mj - maaa at Joint t 

Xj(t n ) - displacement of maaa 1 In the X direction at time t n 
£-* R Xl^ c n* " of the resistances at Joint 1 tn the X direction at time t n 

The twenty simultaneous equations will be numerically Integrated by use of terms through the 
Tlrat backward difference of the difference equation: 

MW ■ *MM - MW* <M>*[i + ev+k* + ^ + ..]*>«!> 

The results will be checked by use of the fourth order Kutta-Glll Integration procedure from 
the Whirlwind subroutine library. 

We are planning on running solutions for a large number of loading conditions as 
well as several conditions of support. 


Charles W. Johnson 

ANALYSIS OP A TWO-STORY STEEL FRAME BUILDINO UNDER DYNAMIC LOADING 


With the program discussed in the last quarterly report, we have oarrled out a large 
number of calculations for a wide range of resistance functions and loading conditions. We 
a loo decided to make some small modlfloatlona In the origins, program bo that we could Investi¬ 
gate anothe: group of resistance functions. This new program haa now been written, and we 
have run solutions for most of the resistances. 


Charles W. Johnson 


RESPONSE OF A FIVE-STORY FRAME BUILDING TO DYNAMIC LOADINQ 


During the last quarterly period, wo have run a number of analyeeo or tho building 
frame using loads of varying Intensity and duration. We have also checked one Whirlwind aolu- 
1 J°n " uh » de »“ calculator computation. The magnitudes of the amplitude of these two aolu- 
tlona agree very eloeely end there la only a alight phase shirt between the poake. 




, . * description of thla problem la glvon in the proceeding reporta fil along with e 

PSl ualn* 0 a"dirraMLI n 2iuiJ? n ' 0 r Utl0 ? outllne of * numerical approach to the aolutlon 

(2J using a difference equation formulation for variable spacing. Since then, a number of 

EKnEi^-XaSSS: hir XT “S' >m . n Che C,rd Calculator end some 

definite conclusions have beon made. Alao, the question or the stability conditions (the 
behavior 01 an error propogatlor.) haa been investigated and will be reported on here. 

The general problem la 


&[*<»&] 


( 1 ) 

In OWxdL, and Initial and boundary oondltlons 


K ( *)> o 


The difference 

equation 

(2) 

/He j 

m 

where 

Si ■ 

j- 


1 

X* 

O 

J- 



, o<t 

- * 

, *** 

XC*,0) - O 

t 0 ■* L 


n 'j 


T*Kj «V- 


WT* 

k(ik. 






i^~r > - J - 


x:) 


and 




f_) 

K„ = k-fx-.) 

The boundary oondltlons are now 

(5) 


x« 






■°'j > *j'j j j> 


U -ij- ' V !,~ * bQ/*. 


Slnce the solution of (2), with variable coefficients, le practically Impossible, 
the ease for constant coefficients will be discussed. It Is easily seen that a surriolent 
condition for ( 2 ) together with ( 3 ) to be stable la that Its domination coefriolent formula¬ 
tion alao be stable [4J- Thus, stability conditions of Che problem 


(A) 


Hr m • J 


* 3V„ - * 


V,,-. , U,.,,-- 

will be discussed.’ The property A+B+C-l la to still hold. 

The term stability ao Is used here refers to the manner In which an error propagatce. 
When an error or the type studied here occurs, it would be due to round-off or gross mistakes 
and consequently can be represented as U fi -C/ s . [5]. the kroneekur delta. This would com¬ 
prise the Initial condition of (A). The Afror anTa: Satisfy (A) together with the boundary 
conditions whloh would be homogeneous as far aa the error la concerned. Further, the error 
Itself Will have a product form of aolutlon. Thua, suppose U n m • N n M m , so. 

(5) -/Vn - Jg.rtH.tV ’ t Afty. = A 

The aolutlon of the first equality la M m •■ It la Immediately apparent that If 


|A| > 1, IAI . 1, IA I < i 
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the error *111 grow exponentially. linearly, or decay exponentially, respectively. The second 
equality In (5) Is 


N: 


a.i 


provided C * 0. Setting 

£ - j + z/J ' 


£ n, e £ V. - 


( 6 ) 


The solution [5] haa the form 

N„ * ( £ ) * [ a «->w-^ e 

where a and b are constants^ Complex va .-a of ^ are not excluded Applying the boundary con¬ 
ditions and considerable algebraic manipulation, may be determined from the relation 

(7) (Si's * N )/"«-* *« ~ £(*/** i*)] ’ * 

where d 2 - 4. Thus. << for r . 1, 2. 3. N-l (the value r - 0 yields s zero eigen- 

function). l A;«o, the case 

Cot l* M « ^ J f j't) 

la another perfectly good zero of (7) and when SuLdt!tuted Into (6) glvea/\» 1. A total of 
N eigenvalues sre found as would be expected. 5* It.* /Al ? | <•*- 

(0) | If e Jl^Ac | 4 I f- ->// f. 

It is Indeed sufficient ror purposes or stability to require tnet ell coerriclenta be positive 
or zero and thus Less tnsn unity: 

|15 ♦ a | tf. i v^‘ ( < IVffltd’l. 

It will be noticed that for some va -ea of b<0 and A * C (8) will still hold for all e^,. 
Likewise, some negative vs ues of either A or C ere not excluded, 

A second and even more difficult question to answer Is tna> of convergence. Under 
what conditions does the difference equs'Son sc.,-lion converge to the solution of the differen¬ 
tial equatlun" A complete answer to this question Is unknown. 

A Numerical Example 

A numerical example has been worked out and the results will be described here. Let 
K - x + b, L - 093. b ■ 2. Q • . In problem (l). Also, let us adopt the rule that all coef¬ 
ficients In (Z) are poattlve or zero. It will be seen mat spacing In the x direction may bo 
taken according to 


will satisfy the above stability conditions. Table I gives the manner in whloh the spacing 
was chosen and the resulting coefficients. Constant spacing was resorted to Tor n>7. As s 
direct result or this, the J values had to te decreased In order to maintain posltlvenesa of 
the coefficients. 

The next step Involves combining tne difference equations ("telescoping”) In such 
a way as to obtain s constant J *»i o Tor the whole Interval (0. 093). The method for doing 
this la described 'n (2). The resulting difference equations are of the Torm 

lA-el* * C A e,. U , „ * A, y n ‘Oj N- 

& 

Those new coefficients are listed In [3] Tor nd6. It will bo noted that (9) reduces to (Z) 

Cor n • 7. 8 and 9. Also. s„ - 0 for n?7. U, n u,, ,. U,, are all readily 

calculated from Table 1. ThQa: lO.mfbV U.*f6». u 13."*-04 ,re 811 re80 “ 5 ' 

" ( * it ij. A ir jt 



} 

* m 

) 


) 


) U *. - 


)V„,~ 
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and similarly for U.. , 

are used when appropriate: 


U 


lZ.m+64. 


U l4,»»64' 0f cour8e U» boundary conditions (3) 


If one were to program this problem at this stage and run it off on a high-speed 
computer, then compare the steady state with the analytic steady state, on-- wo d find trial 
the Tormulatlon as It has been set up here so far does not converge to me analytic steady 
state solution [3J. It la believed that the reason for this 18 that tne left-hand boundary 
condition (J) which reeds in the heat la a very had approximation to the actua. boundary 

!s n ? e the “'"Eblarlty at x - - b - - Z la so near A way or getting around this 
difficulty [ Zj la to require the condition 

no) £ u( 

That la to eay. bv Integrating over the Interval (O.L) numerically, then adjusting the boundary 
value so that (10) Is approximately satisfied, this dirrerence equation formulation may be 
made to converge to the proper steady state solution. An alternative method la to adjust the 
Q so aa to obtain the proper head In-flow (3 J - in a very similar problem [")] variable spacing 

across the entire interval (O.L) wat used and again there was dirriculty with convergence 

even though ( 10 ) was enforced. It la believed that this waa due to the raot that In spite of 
the higher derivatives present In tne truncation error associated with (Z) being very small, 
the exponentially Increasing hj tended to make them significant. 

Table II presents some or the Intermediate and Tlnal results of the example problem 
used here. The teat Tor the dominance of the steady atate la very simple. It will be remem¬ 
bered that in the analytic solution there la a Qt/1. term. If one applies a forward differ¬ 

ence operator with respect to. t (with spacing constant Jk) to the steady state analytic solu¬ 
tion (11), there remains the term 


got. 

gd(mjk) - ^ - B 55 - .0178771. 


IT one applies this operator to the CPC va.ues of U , one should expect to get very 
this numerical value over the entire Interval (0,L)v""Thla teat waa used to determine 
steady atate dominance In this example. At m - 4480 

nearly 

the 


AU 0,4480 " 

.01902. 

AU N.440O * - 01 7“- 


whereas 

AUq.0000 - 

.01790. 

flU N.8000 ' • 0, 786 


and 

AU 0,9152 • 

.01788. 

flU N.913Z * 01787 ■ 



The analytic steady atate values 13J- were calculated from 

(11, V - W ^ - i 

* l ~ * f: ’ 

The per cent error for m - 9132 ranges from -Z per cent up to + 15.3 per cent. Near the 
singularity the values are worst S3 would be expected. A question '-hat one might ask at this 
point: la this dlacrepency near the boundary due to the proximity of the singularity or to 
the variable spacing or both? Two separate problems might be posed to answer this. One. 
with no singularity but with variable spacing; second, a singularity present but the spacing 
aa near constant aa stability conditions will allow. 

It la felt that this method of combining the difference equations ("telescoping") 
holds considerable promise for variable coefficient parabolic partial differential equations 
since It can combine the higher accuracy of a fine mesh with the rapid progress of more 
coarse meahea once the coefficients have been calculated. Considerable experimenting was 
done with variable spacing and It aeema that one should avoid using such spaolng tr possible. 
Various rorma of the boundary conditions were tried other than a difference formulation and 
found to give much better results even with a singularity In the neighborhood. 

The author la very grateful for the saatatance of both Professor P B Hildebrand 
and Dr. P. M . Verzuh Tor their aaatatance. advice and encouragement with this problem. 

Thanks should also go to the Statistical Services Starr of MIT. ror tholr cooperation In 
the use of the CPC. 


Philip L. Phipps 
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Table II 


Difference 


Table I 

Equation Coefficients 


Comparison or Analytic Steady 
(m - 4980, t - 1120 ) 


1 1 

2 3 

3 7 

* 15 

5 31 

6 63 

V 127 

8 255 

9 383 

10 511 

11 639 

12 767 

13 895 


calculated 

U 

n,m 


6.43666 

5.91561 

5.2*955 

4.41031 

3.46331 

2.49642 

1.75301 

1.2*094 

.99078 

.85256 

•77625 

.73921 

.72927 


values 

f lven by 
11 ) 

"57BS239 
5.47713 
4.96740 
4.38277 
3.75430 
3.10740 
2.46376 
1.B48J1 
1.30052 
1.03847 
.89382 
.81358 
.77*13 
.76284 


State with i 

(m - 8000, 
- 

calculated 


7.43664 

6.91549 

6.24919 

5.40941 

4.46130 

3.49230 

2.74515 

2.22718 

1.97276 

1.83166 

1.75355 

1.71556 

1.70537 


Calculated 
t - 20001 
TSTueo - 

f lven by 

11 ) 

B.fle'.tv* 
b.46039 
5.9506' 
5.36601 
4.73754 
4.09064 
3.44700 

2.83155 

2.28377 

2.02170 

1.87765 

1.75355 

1.75737 

1.74608 


Heady State 

(m . 9152, 
- 

calculated 

Vm 

- r. i"uc— 

7.75862 
7.23747 
6.57117 
5.73137 
4.78324 
3.81420 
3.06699 
2.5“892 
2.29444 
2.15327 
2.07513 
2.03712 
2.02692 


t - 22881 


f iver, uy 
11) 


6.78216 

6.27243 

6.68780 
5.05933 
4.41243 
3.76870 
3.15334 
2.60555 
2.3«350 
2.19884 
2.11861 
2.07916 

2.O6766 


EVALUATION OF THE EXPLICIT DIFFERENCE PORMULA FOR A PARABOLIC DIFFERENTIAL EQUATION 

The Introduction to the problem and a description of the procedures for optimiza¬ 
tion of the transient response have been given forth In Quarterly Progress Report No. 15, 14, 
and 15 . 

The optimum values of r, and r for the 7 = 1 case (case 1 ) and the - - 1 case 
(case II) were shewn to vary arolmd R = A, depending on the number of grid tflvlslons, M. 

But ror the t - Q J- * 1 case (case III), 'the optimum values of r, and r were found to be 
significantly lower than ?- when Q Is finite. The envelope R of the optimum values, when 
M-»oo, for this case la shown In Figure 1 . 

The corner values 7 . 0 and 7^-., oeflned in the earlier reports may be viewed as an 
Initial Impulse adjustment of magnitude a such thst for case I 


for caBe II 7 1>0 - jj^A + f(0)J 


and for case III 1 _ x Q - jg- [A + f( 0 )] 

where f( 0 ) Is the boundary lfiput function at t - 0 which la 1 In the above casea. The Initial 
Impulse adjustment A Is Introduced to optimize the transient part of the responae. The envelope 
of A., when M -* — , Is shown In Figure 2 . It Is Interesting to note that the A, value approaches 
that b for case I when Q -» 0 while on the other hand approaches that for case II when Q . 

Numerical experiments, aided by the Whirlwind Computer, show that the use of r value 
aa defined by R In Plgure 1 with the proper Initial Impulse adjustment lmporvea the transient 
responae considerably when « Is finite. Improvement Is optimum when r, and r are differen¬ 
tiated deacrlbed by the present procedures. For oaaes I and II the present procedures, how¬ 
ever, yield no appreciable Improvement over the use of r, * r • R - r. For empirical guide 
on the choice of r and A values for casos under similar Koundary conditions, Figures 1 and 2 
are recommended. 

Boundary Input function or sin oit, Instead of a step, was also tried on the present 
procedures. Numerical reaulta showed that the steady state error was dominant over the transient 
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error and Improvement on the pans lain, reiponae nod only minor effect on ihe overall response. 
It aaemo that tie lower tno r value -lioaen, the smaller would be the ateady state error. But 
on the transient part of the response, this advantage may aoon be somewhat offset by an In¬ 
crease In transient error. 

For details of this wrk one may refer to a complete report with Professor S. H. 
Crandall of the Mechanical Engineering Department. 

Andrew T. Ling 
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STABILITY OF THIN, SHALLOW ELASTIC SHELLS 

bololda 1 a^lTrS™ 1 dU ' rerentlal « 1 ‘'atlona of shallow shell •h.ory Tor . hyp. nolle par., 

m nlH&. - fe)‘- * f 

pv'vr •-e-s-jgj.gn*.. ujgv £ 

where w is the deflection In the z-dlrcoUon a..d K is tt •• ,tress i n. ■ ‘ be will cnatdar 

these equations under the following condtUuna 

(a) The shell Is uniformly loaded by p . 

(b) The edged of the ahell at * • '), a a .<1 , - m. n are assumed to have moment 
free support. 


(c) The edge stiffeners of the shell .ire assumed to be rigid :n the direct 
their axes and to have neg.lglble vending resistance in planes tangent 
she 11 . 


Ion of 
to the 


Using these condition!!, a solution of the Equations (l) :» 

(2) w - 0 F . - p o |£*y 

To teat the stability of this solution, we ass..me new solutions of the form 
(>) w-0+w,F-.p p ||xy . # 

substitute these Into (1) and linearise the differential equotlonn fl] . We then will con¬ 
sider two cases of the stability of the ahell: 

(A) We consider the edges or the snail to be "loe.ued" in the displacements given 

by the linear membrane solution ( 2 ) and then we attempt to rind the magnitude of the load, p , 
which cauaeti buckling. That la, from (2) we may determine not only iho deflection, w, but ° 

also the displacements, u and v, In the x and y directions by moans of the stress-strain 

relations. Therefore, add to the conditions (a) - (cl the condition that, on the edges, u 
and v must be equal to the displacement!' derived frum (2). 

(B) We put no further restrictions on the shell. That la, In contrast to A . we 
allow u and v along the edges to vary from the linear membrane state, compatible with (c). 

The boundary conditions corresponding to Cose A are 

(4) Xr •>* •' *tS ‘ * ° 

^ r <5 ll U/ - f‘£> , 4 > xwa Qyf CJ cu 

where analogously to (}) u and v are the deviations from the linear membrane state. The 

boundary conditions corresponding to Case B ore precisely similar t" those of Case A except 
that the restrictions on u and v are absent. 

To find the buckling load, we assume series expansions for w and $ , satisfying the 
boundary conditions, as follows'. 


S>- 


$ = IT tT' 0 *-- r - f'- -p 1 + 4'y 

Now in Case A the restrictions on the eoge displacements effectively restrict the first deriva¬ 
tives of the stress function, J , through the streaii-stra In relations. For this reason we 
must have the B nn xy In i(j) Tor this term, still enables us to satisfy the explicit boundary 
conditions on J In (h 1 and at the name time maxes lr. possible to satisfy the conditions on 
u and v. However, in Case B J> and Its first derivatives are completely arbitrary on the edges 
30 thot we may without loss or generality take Bqq • 0 In this case, 

The boundary conditions nr Case A make r convenient to express tne problem In 
variational formUj. Having done this, we substitute (cl Into the variational expression, 
perform bhe Integrations and tame the variations and so end up with a doubly Infinite plus 
one set oT linear algebraic equotlonn for the A^. B mn and B 0[) . 

In Case B the lack of restrictions on u arid v make It Inconvenient to use varla- 
tlcnal techniques. Instead we aubst'iut" ( r O into the linearised dlffereitlal qustlona for 
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w and J and than apply rinlta Fourier Transform technique*. OoEn« thill. 
Infinite aet of linear algebraic equations In me A and S ,. Theae are 
the equation* for Caae A except for those teres depending oHTigQ. 


«e get a doubly 
Identical with 


In 


either Caae A or B now nave an Infinite characteristic T* 1 *** 


solve for the characteristic roots p . The smallest root gives u« taucklln* p^. 

In order to solve the problem. He consider the aeries (5) terminated " 1 “*“ 

solve the resulting finite characteristic value problem, be then Investigate *** j’° n *****nea 
of the solution as He add more and more terms to the approximation. The form of the algebraic 
equations He have to deal Hlth la 

(6) 


3x - 




where S and S are symmetric matrices and S la poaltlve definite. To solve this ayatem Hhen 
H and 3 are finite, «e first reduce S to Its canonical form and then diagonalize H. This 
amounts to the simultaneous reduction or tno quadratic roras to "uma of squares which may 
alHu:'B be done If one of the matrices la positive definite ( 3J. Actually, the Equation (b) 
depends also on tHo parameters not explicitly shonn: 


(a) a/b the ratio of the linear dimensions 

(b) h/c the ratio of shell thlckneaa to ttiell rise 

Computations are non under nay on the WhlrlHlnd I computer for various values or 
the parameterB. Thus Tar no more than tnelve terms have been used In the approximations and 
Tor values of h/c Hhlch are not too small this number or terms appears to give quite good 
convergence. 


Anthony Ralston 


References: 

fl] Machine Methods of Computation and Numerical Analysis, Quarterly Progress Report No. 15, 

P. 12 , March (1955). 

[ 2 ] E. Relaaner, Jour. Math. Phys. 32, Nos. 2 , 5 (1953). 

M A. Kecslsr, quarterly Progress Report, Solid State and Molecular Theory Group, p. 15 
Ootobcr (195*5. 

STEEPEST DESCENT ANALYSIS OP RESISTIVITY DATA 


It Has round that many resistivity analyaes attempted using the Neuton Method gave 
unsatlaractoy results tlj. The error "aurracea" Along Hhlch the Iterations carry the solu¬ 
tions are generally non-llncar enough to cane difficulty, usually extreme overshoot. It 
hub necessary to compensate for thla. during the rirat dozen or ao Iterations, by reducing 
the mremoter changes made to about five per cent of that predicted. Only when very near 
iheaolttion (fits better than one per cent at all Kernel points) Has It occasionally possible 
to use the entire prediction. Another unaatlsfactory feature or the Kenton solution Has tnc 
necessity or solving a matrix to rind the changes at each Iteration. The matrix solution 
could Introduce lnaocuraoy in a variety of waya (two variables nearly equal, giving degeneracy, 
or one variable orders of magnitude larger than another, or much closer to Its solution). 

Theae conditions result In Indeterminacy and Hlld predictions, or In singularity Hlth the 
consequent machine difficulty. 


Apparently the Newton method was not too well suited to the Rinds of functions In¬ 
volved. and It decided to Investigate the deepest descent solution for sets 0 f non-linear 
algebraic equations. Thla la designed to use more data points than variable*, giving a lesst- 
aqusre fit. It does not require solving a matrix. In addition, a way .as round to compen- 
sate to a certain extent for the large dirrerencea In magnitude of the variables, a result of 
their dimensional disparity ( f - resistivity, d . length). Thla procedure was found superior 
to the Newton procedure and Id described below. v 


Reconciling the notation of the previous report with that of Householder: 




fj(h)- F(h) * f,i 

i* ’ ttVi " i 3 

r, = -T^r- - 7 



t ' 71+! - 

*/(?/)» X [° ( (t-rt.lr f S ,'Panama tj 

Proceeding with Householder's notation 


X, - X 0 - X « 

and. In general 

x st> ■ X, - U. 

( Ue,* 
u, r <.*</)' 
th (u,).. (Vfj, 

the 1 components of v* after a lteratlona, and X la auch as to minimize 9 . In particular, 
for a three-layer caae, f 1 and known, /j s f 2 . d 2 . - dj 

VJJ * (// f (**/i y.r 3 '* v 

V4>,. (**/>£)/[ W/apy*. 1**44)'’} * 

Aa mentioned previously, the dimensional disparity of the unknowns results In aome 
of the u, being much larger than others. Thla, In turn, causes the variables with smaller 
derivatives to be relatively Ignored In the iterations until those with larger derivatives 
are very nearly at a minimum. Such la the caae when the second layer la thin and very reala- 
tlve--the resistivity la virtually Ignored until the end of the solution. Compensation can 
be achieved by using. Instead of the Uj given above, 

u t . —Llf_ 

This worked satisfactorily. 1 • 


The optimum value of X can-be.round ajactly only by trial end error. For any 
function, 0 . It can be approximated by wp, or i™, etc., each aucceaalve approximation being 
exact for a higher order power series (Newton Approximation). However, If the function la 
not expressible In a power series expansion, then the approximations are not necessarily valid, 
and aucceaalve Newton approximations may become poorer. Such, for Instance, la the case when 
non-integral powers appea: In the function. 


Examination of the error surfaces showed that negative second derivatives were not 
uncommon when at some distance from the minimum. Near the minimum, however, curvature la 
always upward. The calculation of A was handled as follows. From the first estimates. 


s<Ji' 5 


V- 


the function «> was calculated, aa were 


and 




the first and second derivative* of * along the gradient path. Then new trial values, 8 , of 
were calculated from , 

i t - * - c, (+•/*;) <x. 
x. - (*;/?: )u. 


Obviously 2y is exact for a linear function, and 24 for a parabolic runotlon. There la a 
dlfference u between theae two approximations that la of considerable importance when the 
error surface, as here, doea not go to zero because of errors in the data. It lo seen 
that a linear extrapolation used under theae conditions will always overshoot. 
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such mat 


Three-Layer Analysis of Three-Layer Kernels 


o, * °0,)- 

A value or c was then round uy In erpol.. '.on w< ' • *‘>uld '.nlmlze a. Because of me behavior 
of the error curves, the value c could • tta! .»t 1< m dare manner when near th solu¬ 
tion as when at some distance way. tt fe uij.i : elnte^ol stt on^to^ft 0 
at the beginning of the eel - , end S psret - ' I lallort to 0 • mlnlous. at the end. 


That la. 


tglnnlr.g 
near • h 


beginning, assure 

c r 

2 . 


* 


£>. - A 


at the end. as.mme 


-I 


a ,u uae 


V 

O -- 


J' 

J 


changing fro. > linear 


c C' 

i*C 

- JL 

parabolic when 


n, \ M, - 4>.) > o/l <P, n>,J V C [tj - d>J 

6 ! (-&) * gjj (a-Cf) C- t ) 


In tnc cases analyzed, values of o have general been betaken 1 2 aid 5 2. 
Proceeding aim the Iteration men. 


V*- 


etc. 


an “ir.terpolat ton within an interpolation" 
wo "du j " iterations for every true ltera 


When this worn was begun. It was felt blur auen 

procedure might he lnefriclent because of the two "dummy" Iterations for every true Iteration. 
In practice this lias not proved the case, since er: ro dlmlnlah a*, four to five times the rate 
of the ordinary steepest descent procedure on the same data. For more regular error surfaces, 
of course, It la protablo mat only on- rind of Interpolation (linear or parabolic) would ne 
satIsfactory. 


Results 

Most cases analyzed here were t.,oee In wntch difficulty was encountered In the 
Newton analysis, (Vozori 1> ). Ir addition, two other "difficult" cases and two of the 

field cases were re-analyzed. Comparison or coses , ii. 10. and li with the results of the 
Newton analysis Shews t.ne greater power of the Steepest Descent analysis to reduce the error 
of the solution. However. It will it noted that the solutions It. me two rield canes, B and 
C, arc not, significantly different fro; me Newton solutions as might he expected. The thick 
second layers ma.ee then "easy” coses. 

Cases 10s aid 10b anow the results or some of the more common Steepest Descent pro¬ 
cedures or the data or case 10. Both is- the nor-norma 11 zed u, . Case 10a was done changing 
one parameter at. a time. Cose I0t. waa solved changing all parameters simultaneously, differ¬ 
ing from the rest or the analyses only lr. the non-normallzed u,. Both methods are slower 
than the normalized gradient group descent. 1 


Conclusions 

The normalized gradient group descent analysis presented here possesses two advan¬ 
tages over other methods .Bed for the resistivity problem. First, oy virtue or Its Steepest 
Descent nature, it has a better chance of converging to a solution man does the Newton method. 
Second, the norma 1 IZBtlon and type of Interpolation used allow mote rapid convergence to a 
solution than do the standard Steepest Descent techniques. 

This technique probably should be applied to cases of more than three layers. 

ft is again a pleasure to acRnowledge the advice of T. R. Madden. 


Reference: Keeva Vozoff 

[l] Machine Methods or Computation and Numerical Analysis, quarterly Progress Report No. IS, 
p. 55, March (1955). 
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Three-le.yer Analysis 


First Estimate 


ol' Pleld Kernels 
Sol At Ion 


r : r 2 

<*t d 2 

1 . .0 

,H A. 

1 . . 2*1 

.3 U. 


r, A A A 


.2 1. .025 .2 

I.aIO" 11 .1V« 2.72 1.3X10* 5 

.2101 1. .251 .2101 

3.3xlO~“ .272 1..X10'" 


5.1 Institute Courses snd Seminars 
COURSES ON DIGITAL COMPUTER CODINO AND LOOIC 


Course 6.555, Introduction to Digital Computer lodlnv and Louie , a discussion 
of selected topics In programming, logical design and applications of large scale digital 
computers, was orrered by D. Arden at MIT during the sprite or 1 • . ! e course included 

tne solution of a programming problem on a simplified single address computer simulated oy 
Whirlwind I. Among tne probleme solved oy Chit clase members were the solution or simultan¬ 
eous linear equations, Integration of differential equations, and tne economisation of 
power aerlea. Tne total enrollment was 55 seniors snd graduate studenie (from 1,0th tne 
engineering and Industrial management curricula). 

The Digital Computer laboratory Comprehensive System (CS 11) programming course 
was given once during this quarter. The couroe Includes tne following topics: relative 
addresses, temporary storage, rioatlng addresses, preset parameters, programmed arithmetic, 
cycle counters, burrer storage, automatic output, poat nortecs, ami 1 .lilpsas conversion. 

The text for tne course la a programmer^ manual written y atarf n- i.ern of tne S and EC 
Oroup. The 26 students enrolled during this quarter represented tne following groups: 
Department of Business and Engineering Administration, School of Industrial Management, 
Department of Nuclear Engineering, Spectroscopy laboratory, laboratory for Nuclear Science, 
Naval Supersonic laboratory. Chemical Engineering Department, Physics Department, Solid 
State and Molecular Theory Oroup, Aeronautical Engineering Department and Lincoln Laboratory. 

NUMERICAL ANALYSIS 

This year, for the first time, Professor Htldeora.d gave a second semester course 
in Numerical Analysis, MA12. Tula semester tne context or the course was made up of the 
following topics: Least-squares approximation, smoothing of data, quadrature rorimilaa or 
Gauss and Chebyahev types, harmonic analysis, exponential and trigonometric approximation, 
determination of periodicities, Cheoyshev approximation, continued fraction expanalona and 
rational-function approximation, numerical solution or partial differential equatlone. 

MACHINE METHODS OF COMPUTATION AND NUMERICAL ANALYSIS 


During the last quarter the biweekly seminars for the project reported In Port I 
were completed and brought to a close for the academic year. Tne nigh points Tor the 
quarter were two lectures presented by guests from the National Bureau of Standards and two 
lectures Jointly presented by the Meteorology Department and the Mathematics Department of 
MIT. 

The guest lecturers were Pr. Milton ADramowlt* and Dr. Frans Alt. Dr. Abramowltz 
spoke on "Coulomb Wave Functions”, a topic which Is Important to tne phyolclata of tne 
project, as one can tell by referring to the comments on Oroup Activities prepared for the 
previous progress report. Dr. Alt spoke about the type of problems which na yet present 
overwhelming difficulties (storage demands, etc.) for the modern electronic computer. 

The Joint lectures by the Meteorology and Mathematics Departments brought results 
which are summarized In the report entitled: “Multiple Prediction" ( Part 1, Section 2.2 
above). These lectures presented the multiple prediction techniques currently being applied 
by the MIT Meteorology Department as well as the Matrix Factorization work of Professor 
Wiener which should lead to s scalar solution of the multiple prediction problem. These 
two points of view (cf. problem 155. Part II, Section 2.2 and reference (5] of Multiple 
Prediction) revealed the need and potential availability of a scalar solution to weather 
prediction In the multiple series case that has previously been available only for single 
series. It seemed profitable to look Into the matter further and this waa done during some 
meetings of tne probability seminar of the mat.nematics department, the numerical analysis 
discussion group of the project (see Section 3.2 below) and during further conversations with 
Professor Wiener. 

3.2 Group Activities 

NUMERICAL ANALYSIS DISCUSSION GROUP 


The usual membership and purpose of this group nave been recorded In previous 
reports. The majority of the discussions during the last quarter were held for the purpose 
of obtaining the multiple prediction solution (see Part 1, Section 2.2 anti Section 3.1 above) 
snd consequently were attended mainly by Professor A rtnsnd Siegel, Physics Department Boston 
University, Dr. Joseph 0. Bryan. PIC Staff at MIT, Mr. Peter Hanna and Mr. Robert Miller, 
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Meteorology Department, MIT and Dr. Bayard Rankin, Matrrmntlca Department, MIT. 

COULOMB WAVE FUNCTIONS 

The programming or Coulomb Wave Functions, Doing carried out Jointly by some 
phyalca department members of the Committee on Machine Methods of Computation, la still 
In progresa. 

NUMERICAL ANALYSIS LABORATORY 

The nupervlnlori or the Numerical Analysis La moratory, which la open approximately 
0 hours a week In conjunction with. Professor Hliaeorand'a courae W* 12, and the grunlng oi 
homework problems has been done by M. Douglas Mellroy, Philip M, Phipps,and Anthony Halston, 


PART II 

Project Whirlwind 
1. REVIEW AND PROBLEM INDEX 


This report covers the speciric period of March El, 1955 to June 15. 1959. 

During thla time 74 problems made uae of 3?2.8 houra of the 546.5 houra or Whirlwind I 
computer tlmi allocated to the Scientific and Engineering Computation (S and EC) Group. 

The remaining 813,7 houra or the allocated time were uoed for terminal equipment teatlng 
and calibration, demonatratlona, tape oonverslona rot Lincoln laboratory and various 
lnter-run operations not logged to specific problems. 

These problems cover some 18 different nelda or appllcotlona. The results or 
E7 of the problems have been or will be Included In aoademlc theses. Or those, 19 repre¬ 
sent doctorate theaea, 3 engineering, 7 master's, and one bachelor's. Thirty-three of the 
problems have originated from research pro.lecta aponaored at MIT by the Office or Naval 
Reaearch. 


Two tablea are provided aa an Index to the problema ror which progreaa reports 
have been submitted. The rirat table arranges the problems aooordlng to the rield or 
application Indicating the source of each problom and the per cent of the WWI machine time 
consumed. The second table attempts to arrange the reports aooordlng to the principal 
mathematical problem Involved in each. In each table letters have been added to the prob¬ 
lem number to Indicate whether the problem la for academic credit and whether the problem 
Is aponaored. 

It may bo seen from table 2-1 that problema originating from ONR sponsored research 
used 60.ol per cent of the total WWI problem time. Thla 60.61 per cent Includes 11.78 per 
cent used by the Digital Computer laboratory staff, 9.67 per cent used by the Solid State 
and Molecular Theory Group, nnd El.30 per cent by the Machine Methods of Computation Group. 

A Data Reduction Symposium was sponsored at the Digital Computer laboratory by 
the Servomeohanlamo laboratory. A report on thla symposium la Included under Problem 1E6. 

A routine that computes and displays some 15,000 flve-dlglt numbers on 145 rramea 
In about four and one-half minutes was uaod In Problem 224. 

The Whirlwind programmed-marglnal checking facilities have beon expanded to Include 
additional terminal equipment. 

Reliability flgurea are provided for the Whirlwind computer system covering the 
3>-week period beginning on 28 September 1954. During this period, the average uninterrupted 
operating time between failure Incidents was 10.8 houra. 
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•- -.:vr CODING AND APPLICATIONS 


2.1 Introduetl.- 

Progres: -* -• ’ ;.;-'.tted by the various programmers are presented In numerical 

rder In >.t « » >-—ary report presents the combined efforts of DIC Projects 

O end 01- . -T: undertaken by members of the Machine Methods of Computa¬ 
tion Group have .<• - ,>.-::on 2.2 or Part II to avoid duplication or Part I. 

Suitable cross -.luded In Section 2.2 of Part II for completeness. 


277 . 278 , 280 - 28 *. . 
Tlrst time. Eleven -. 
have beer, complete:. 


red :elow, thirty-one (22o, 256, 259, 201.-206, 207-275, 
rr:resent new problems that are being described for the 
' ■ . 2- , 201, 268, 282, 283, 286, 292 and 29*) 


• ai .es —-it.- - - ee- ae- up to provide the reader with a convenient Index 

to various ln-eres- , • • •.» pr.'blema. Table 2-1 lists the problems accord IngTo— 

i" #, r :s •' "■■■' -a-.es tr.e source or each problem and the percent of the 

ar.fl r.c Group’s .. *- "2.5 r. ..-a) used uy each. The mathematical problems and 

procedures ■ . -• . -srl.-us rble-a «r- tabulated separately In Table 2-11. 

■ oott ables. pr > • .reflx-u asterisks (a) represent work being performed 

By members of the Ms: : e Me”..: a ' f Co-putatlon Group. 

- . . l^tlcra have ee ,:lti . tne problem numbers to indicate whether the problem la 

Tor academic cred.t aruj . ---er - s sror,sored. The lettern have the following significance: 

A Implies -r.e ;r.:Ie- :a NOT for academic credit. Is UNsponsored. 

B Implies tne prctlc- ^ academic cred'.:. Is UNsponsored. 

C Implies tne problem •» NOT for academic credit, IS sponsored. 

D Implies the problem a for academic credit, IS sponsored. 

K Implies tne problem la sponsored by the OKR. 

L Implies the problem la sponsored by Lincoln Laboratory. 

The absence of a letter Indicates that tr.e problem originated within the S and EC Group. 

0n Sy8t *” 8 Engineering has boen Included In the Appendix of this Summarv 

2.2 Problems Being lolved 

100. COMPREHENSIVE SYSTEM OF SERVICE ROUTINES 

The Director Tape Program 

several 0rd ' r t0 ™ 

oaaca the 1,181 *" 8 — of 

In particular, a programmer may .U.r th2 Hn^ control TlLi 2J r ! e J? r COpe 1,111 be ““ d - 
computer la not stopped artar read-ln but Instead bilglia JhS comDuta^ ry 8P * 90 th8t Ch * 
may terminate hla program by reading in the next the computation; or a programmer 

If a director tap. l.'not being used, the progf-amS^cJn I* 0 "*"* the <=°'" puter - 

reatart button), rrom the performance request*™ the ?i™t rB (pre8B 

iSa r r S b,^lid n ed b ? tl0n) ' r th * 8888 » • dlr eotor & £*J£l ^^0^. 

performance requests in b^th^r r thi e sbove C csiei*£lll‘’b^the^a^J"* modlried 80 the 

The following new word, are being added to th. director tap. vocabulary, 

(1) al 1 switch on, 

(2) al 1 switch off. 
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The “.r “■ m te to ~ -«.—drum. 

8 Binary tape or a C l occur after 

lerred to thTSiJJcto"^ ^grj^lh^fU^'fo^rSi^tS?"^ ’ , ‘ U ° 8 l ™"" 
ferred imm^la^’to tS'&JE'jSS ™d"5nf 0 "° rr *' ' h '" °°" ppt8p •«*«•» win be trane- 

Recqrd_ln S _ot L P lad Table on th. Bufr.e n.-.- 

18 left of nMt,n * 8dap888 * 8 

can be printed out with £?mboi“'addressee!'"* ° r p08t - m0rt< ”” Programs In which Inatructlona 
106 C. MIT SEISMIC PROJECT 

in "Detection of R^t l^K?^ 

!S*S^S^EnSSK~K2^K‘sSsSr s * s S.»?W*4B SfPSSSs 

of th, MIT^DepartSnr'o^Oeology’and^eophyaleaf" Sln ‘ ,8 °"- »* Treltel, and I. c.Inan 

R e £ erenoe Oeoiog^nd' O^ophyalca 

I. Wadaworth. Roblnaon. Bryan, and Hurley -GEOPHYSICS. Vol. 18 , No. 3. July 195, 

120 B ’ N - SSJS ^ ""*« ° F REJECTION INTO HJOH-TEMPERATURE, 

ponent, called 8n ? "ierofLrmopr«aor ? !^n^hlcr^nerr’se^ 0 P ? t * ntl#1 en0 turbine com- 
gas stream la brought about by the evapora”on of ??n Ti" ,' n 8ta E nBtl on pressure of a hot 
region of the riow. The concepts underline the " 8t ; r Sn J* eted Into s hlgh-velocltv 

ln t f^*'e h c °7P ar8tlv *>y recent won.: }?! the f'eld ” J® rot 'j'rmopressor are an 

gas turbine oycle la sru,logous to thsl of the co£en2e? Tn ^J’n^am'^owerTa"? r “ n ° tl0n 

d?fri tU rb‘^*l"to 1 a* c lrcular*iluet , of 0 varylng*d?ameter"termlnated*r* t " S th8 8 * h8 ' 181 « 8888 ^rom 

ohr^ Ser i,"?i Ch recov8r8 kinetic e«rgy or ^i no'^PS^ d 'a I 80 '; v "’ttonal conical 

p ^ r ®’ . At t»>« entrance or the duct, special lnjeetorl del d J««harglng It to the atmoa- 
are In turn atomized by the rapldly’moving gai .tream? d8l,ver mlnul: ® >»« of water which 

taneoua thenscklynTmlc ind 8 dJ™mT’ , i'fr, 1 |'" Sr^if'^aporaMon 8 ” thrfl 1 M°“ l by th8 ,ln,ul - 

atagnat lon pressure across the device. K^ner descrJptto^ o? ft, 1 ? - ,b “ t 8 rise lr° 
earlier reports, beginning with Summery Report No. 32 ? P^uMh Su^t^rf'i^.™ 1 ' be fo , " d ln 

lntlmat,Iy T connac*ed r wlth r the r 'detarmlnatlon C of*performance P charactariitlca°^f*the P dev 8 oa 
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under nil conditions - s s comprehensive one-dimenalonal analysis of 

the process. This *n«lysis - aneous solution or seven, nonlinear, rirst- 

order d 1 fferentlal equa? *.or.> . 

During the pan »*. lor.a were reaumeri using a new Whirlwind program 

(see preceding nummary Report . “ ese Imputations were divided Into two categories: 

(a) those devoted spe ' .erotnermopresaor performance, and (L) those 

connected will further tee . - - and with studies of the various methods or 

numerical analysis which II ?. 

To Illustrate • '.e •• > « typical calculation falling within 

category (a) above will be : The problem waa to calculate the auperaonlo 

' leal performance of an teretf-er-cpre s r sonalatlng of an 11" diameter shell, seven 
Tern In length and tem-lnn' :: ffuaer having an Included angle of five de¬ 
crees. Within the outer duet k( — »• from Its entrance, an axlally- 

avTrretrlo plug, having a para. . variation from zero at Ho noae to a maximum 

or 4 ,lil-inohes three-and-ane-na'-:' treasi of the hoae, waa to lie along the axla 

of aymmetry of the duct. A normal s-..-.-- fro - auperaonlo to iubsonle flow waa desired 51 - 
hea downatream or the d The sequence of calculations were fully sutomatlc 

and Included provision for vsr> lng tr.e -er* size to weep truncation error within a 

reasonable bound. The Singular solution critical) was automatically determined by an 
Iterative procedure wnlch ad.'ustcd - to I - 1 a I Mach Number until a transition from aubaonlc 

ptraonla now si. consist.- area -co-pllsned. One foot from the duct entrance, 

the calculation of liquid temperature was a :•omu• Ically modified to prevent later oscilla¬ 
tion In tula function. These csl.-1.> isi I the second-order Runge-Kutta method and re¬ 
quired about twenty minutes of computer rime. 

In the oategoi v of further '. sting of ' ne new program, a routine waa written for 
accurately timing the cnloulatIona using me Whirlwind clock. The results are summarized 
below: 


NUMERICA L METHOD 


TIME PER INCREMENT . 


Euler, First-order 0.9b 

Backward differences, first-order 0.98 

Runge-Kutta, second-order 1,83 

Runge-Kutta, fourth-order 5,56 

Forward and Succoanlvc differ¬ 
ences (three-point formulas) 2.68 

Wet-bulb temperature calculation 0.2 


APPRO XIMATE TIME PE R 
WPPE' OPT 
TTT sEcoNds 


r 


8 


These comparative results between the new program and the earlier one provide 
on excellent example of the saving In computer time ths' can be realized with elimination 
o! Interpretive programming whenever possible. The new program la nominally four times an 
rest, and when the wet-bulb calculation routine la required? it la about ali times aU ra" 
(for the fourth-order Runge-Kutta case). It ohould be pointed out that, the significant 
Ml? oalculatlon? Ulb “ lettUtl °" '- 1 ™ '» ■»«■«» to a complete revlalo^"^ 

moat pvor^n;;!^ 

computational time to a minimum while also maintaining adequate accuracy Th«e JtSdlM 
involve a computation of the stream properties over o given Interval us 1 A* a vJrtStv Ar 
increment sizes and then extrapolating to zero Incremen Tlze toA the titict' aS tuot 
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(All quantities 

refer to stagnation pressure rotlo 

• P«/Po»" Bt Z 

- 0.01 feet) 

ACCUMULATED 



TIME OF 

NUMERICAL METHOD 

ERROR 

INCREMENT SIZE 

NO. OP STEPS 

COMPUTATION 


o/o 

roet 

• 

seconds 

Euler 

0.073 

0.00028 

J5.7 

.2 

2nd-ordor Runge-Kutta 

« 

0.00140 

7.15 

13.1 

4th-order Runge-Kutta 

II 

0.00213 

8.23 

29.4 

Euler 

0.26 

0.00081 

12.32 

11.8 

2nd-ordcr Runge-Kutta 

n 

0.00249 

4.01 

7.3 

4th-order Runge-Kutta 

H 

0.00351 

2.8b 

10.2 

Euler 

0.70 

0.00162 

6.1b 

5.9 

2nd-order Runge-Kutta 

" 

0.00371 

2.70 

4.9 

Nth-order Runge-Kutta 

ii 

0.00490 

2.04 

7.'3 

Euler 

1.11 

0.00214 

4.07 

4.5 

End-order Runge-Kutta 

ii 

0.00461 

2.17 

4.0 

Ath-order Runge-Kutta 

M 

0.00570 

1.75 

6.2 


(• The non-integral values appearing here were obtained Trom a curve or truncation 
error vs. Increment size, and do not represent, the actual number of increments uaed In 
computing the ourve.) 


The data above Indicates that the 2nd-order Runge-Kutta method la aIgnlfloantly 
foster than either of the other methods for the ranges In distance and in truncation error 
considered. An unexpected result is that for errors on the order or one percent, the 4th- 
order Runge-Kutta method Is slower than the Euler method. As the error decreases, both 
higher order methods become Increasingly superior to the Euler method. Though It la clearly 
understood that the above results apply to the solution of a particular set of differential 
equations, for a given set of Initial conditions, and Cor a particular Interval, they 
nevertheless Illustrate the Importance or considering higher order methods or calculation 
and, further, of providing alternate methods, when planning a program or numerical analysis. 

It la planned to continue those studies of truncation error In Intervals succes¬ 
sively further downstream. Tills will be done on a limited scale or effort while continu¬ 
ing Aerothermopres3or performance oalcuiatlona. The letter will be carried out using what 
la considered to oe a oonaervatlve Increment size. 

The AerothermoprcB8or development program la being carried out at M.I.T. under 
the sponsorship of the Office of Naval Research and Is directed by Professor Ascher H. 

Shapiro of the Department of Mechanical Engineering. The theoretical aspeots of the prob¬ 
lem treated by Whirlwind I are being carried out by Dr. Bruce D. Govrll. 

B.D.Qavrll 

Mechanical Engineering 
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122 N. COULOMl' WAVE FUNCTIONS 

See Section 2.2 ctP Port I. 

12? P,N. EARTH RES IflTIVlTY INTERPRETATION 
See Section 2.2 oT Par'. T. 


12t D. DATA REDUCTION 

Problem 12o la a very lore? Uata-rednclloh program 1'or joe lr. the Servomechan'sms 
Laboratory. The overall problem In composed or winy component sections which nave been 
developed separately and are now .ring combined Into complete prototype programs. Desorlp- 
tlona of the various component sect Iona I IVS appeared In past piu.-'or ly reports. Al'ter the 
development and testing of the prototypi rams are completes, t«.e programs will 

te re-coded for other, commercially aen:'.able, large scale computera, (probably the ERA 
110?, IBM 701 and IBM 7CO* computers), tor use ny Interested agencies Cor actual data reduction 
at other locutions. The program* are currentl.v being developed ny Douglas T. Ross and David 
F. McAvltm with the assistance of Miss Dorothy A. Hamilton, ilcrvomeehai.tsmo Laboratory atari' 
■embers. Thin work Is sponsored by he Air Ferae Armament I/moratory through DIC Project 
71 Sft. 


The nature of the problem requires no: only extreme auloiietlclty and efriclency 
lii the actual running or the program, tut. also the presence or human operators In the com¬ 
putation loop I'or the purpose or dee la lot. malting and program mod inflation. For thla reason 
extensive use Is made of output oscilloscopes so that the computet- can communicate with the 
human, anti manual Intervention registers so that the human nan comnuntoate with the computer 
In terms of brood Ideas, while 1 lie computer Is running, and have the computer program tranB- 
late theao Ideas into the detailed steps necessary I'or program modiricatlon to conform to 
the human operator's decision. The program which does thla translation and modification la 
■ailed the Manual Intervention frogru.m (MTV). The -os'- rectm version of the prototype 
data-reduetIon program la called the Basic Evaluation Program. 

On June 1 and 2. 1UV. a Data Reduction Symposium was held at the Barts Building 
under the sponsorship of the I'ervomeehanlams Laboratory. The program Included detailed 
considerations of hardware and programming for 'annul Intervention switch Inputs and scope- 
type output a as well ai, the uae or these devices for Incorporating the human operator Into 
the automatic cor.ip . jtlon loop Pom for routine computer operation and Tor large scale data 
reduction. Talks were presented by personnel of the Lincoln, Digital Computer and Servo¬ 
mechanisms I/i bora lories. Demonstrations included the operation of tne Comprehensive System 
Utility Routines, tne programs or Problem 12' and a special Qulded Missile Launch Simulator 
program written under Problem 12n to demonstrate the use or the Ctlaractron tube on the 
Memory Test Computer at Bedford. A tour or the Chaructron and Typotrori manufacturing Facil¬ 
ities al Barts Building completed the program 

The Basic Evaluation program has beer rewritten to eliminate the uae or buffers 
and to Include an nddltlonal act of equations. Some progress has been made In tile rewriting 
or the Manual Intervention Program so that H can ise an arbitrary number of drum groups and 
have facilities Tor enny expanoion. 


The so-called Fart I Progran which has been In use for over a year Is being re¬ 
written for more flexibility and speed. In the process a new (?t,0) square root routine has 
been devised which uses u combination or Whirlwind operations and direct modification and 
entry into the Programmed Arithmetic routines to obtain a square root in •, mllUseconds 
Instead ol the ?1 milliseconds required for the fastest routine In the Subroutine Library. 

i Bo ? 1 ° Input-Translation Program for the ERA 110? Computer Is being written In 
110? code Tor translation by the program or Problem 2-:,. . The routine allows the use or 
mnemonic code, absolute sddresseo, an" Integers to any base, and Includes many checks Tor 
Illegal characters and comb Inst loria. It follows the jor.ventjona or Problem 2ho and Is in¬ 
tended to be used Tor making dorrectlcna to 110? programs while at the site or the 110? 
computer. •* 


1?1. SPECIAL PROBLEMS (STAFF TRATNINO, DEMONSTRATIONS, ETC.) 


D.T.Ross 

Servomechanisms Lab 


As part of the 3 and EC Group activity, 
when necessary to new users of WWI. In the past, 
and discussion with no opportunity for actual use 


a two-weev CS Training Course lo oTTered 
i-hlo course naa consisted solely of lectures 
of the computer. 


»lt.i 
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In an attempt to Increase the effectiveness or the course, WWI time wna msde 
available to students when the course was given during thla quarterly period. The students 
coded a change-making problem, which waa chosen to be simple enough I'or even novice pro¬ 
grammers to complete In a short time. About twelve of the enrolled students van programs 
during the time allotted for this purpose. 

During the past three months, nine groups visited the Laboratory. The affiliations 
or some of the larger groups are given In Appendix ?. 

1?2 D. SUBROUTINES FOR THE NUMERICALLY CONTROLLED MILLINO MACHINE 

During this period, an Electrical Engineer's thesis entitled Automatic Data 
I'reparatlon for Numerically Controlled Machine V o’a waa submitted by ,1 .K.FI uuyot. to the 
Mf'T E1 eclr 1 ca 1 Eng 1 rieer 1 rig Departmenl. This report contains chapters on mclal-cuttlng 
geometry, numerical methods, description of Wwl subroutines written especially Tor the 
milling machine, description of milling machine tape-preparation routines, results or a 
study of computing facilities outside MIT, and suggestions for future extensions or the 
present study. An Appendix contains codes of the library subroutines used during the 
study and detailed descriptions of tlielr operation. Portions or this thesis will be 
available In engineering report form at a later daio. 

Application and further development of the subroutine library wore continued. 

Two subroutines wore successfully tested. One la for more accurate feedrate control 
than was possible with the feedrate-control subroutine previously written and tested. 

The other subroutine la for computing cutter-center offsets given two orthogonal tangents 
to tho surfaoo at the point under consideration. Milling machine tape was prepared for 
rough and finish cuts on a cum for Ctromberg-Carlson Corporation. Another milling machine 
tape was prepared for a spiral to he used In s demonstration of u new numerically controlled 
spin mill at Ulddlngs and Lewis Corporation lr. Fond du Lac, Wisconsin. 

Tile programmers for thla problen were .11.Runyon and T. Nagle of the Bervomechsnlama 

Laboratory, MIT. 


J .11. Runyon 

Servomechanisms Laboratory 


Ml. J AND EC SUBROUTINE STUDY 

An error was found In .lbrary subroutine L3R FU U b. Thla subroutine was designed 
to find nines and cosines of large angles, reducing the angles to flrot quadrant equiva¬ 
lents vy moons of non-Interpreted Instructions, thereby shortening the calculation time. 
However, this reduction was found to be In error In certain cases. The subroutine haB been 
corrected and Is now being tested. 

A seven-point Oauss-lritogratlon subroutine nos been coded and tested. 

Two ro-'tines for the direct or delayed prlnt-out or (?0,0) and (13,0) decimal 
fractions, available In an older subroutine library, nave been brought up to date and are 
now available. 

L3R OS 2, Decimal Integer aoope output routine, has been rewritten to conform 
in the present Hiatus of scope decoders. 


1SI| N. SELF-CONSISTENT MOLECULAR ORBITAL 

Thla project has seen described In previous reports by Dr. A. Meckler, who has 
completed the revisions which he undertook shortly before leaving the M.I.T. Solid State 
and'Molecular Theory Croup. (Eeo uimmary eporl No. "o, p. 1 1). A second revision use Seen 
completed ov Dr. R. Heatet and tested on simple examples. The program la now designed to 
■arry out approximate Hsrtroo-Fock calculations for wave-functions of arbitrary configura¬ 
tion. using the method of symmetry ond equivalence restrictions. In the present program, 
all one-electron Integrals are transformed to the ortnonormal bsoln of self-consistent 
orb' tills after the bast Iteration. The Schmidt orthonorma 1 Izallon process Is usod instead 
of the earlier A technique. Convergence has been made more rapid by an extrapolation 
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method. The program Is being tested In a series of ealouiatlons or gradually Increasing 
complexity, described under Problem 288 . 


R. K. Neabet 
Solid State and 
Molecular Theory 
Group 

155 N. SYNOPTIC CLIMATOLOGY 


Durli.fi the last quarterly period, the objective or the project has been to 
spoolfy and predict , Day Moun Anomalies for both temperature and proclpltotlon using the 
same Statistical methods which the project has been using for the paat two years. Thl3 Is 
o continuation of worts which has neon In progress for the past Tew months. Because of the 
amount of latn used, and the site of the matrices Involved, It became necessary to program 
the arithmetic operations; several programs have been written for this purpose. 


At the present time this particular phase of me pronlem Is nearing completion. 
Further edmmenta concerning this problem cannot be made until all results hove been com¬ 
puted and analysed. 


Programming is being performed by Mlea Elisabeth Kelley under the Eupe.-vlslon of 
Professor Henry G. Houghton of the K.I.T. Meteorology Department. 


E. A. Kelley 
Meteorology 


COEFFICIENT IN A SEMI-INFINITE RECTANGULAR WAVE QUICK 


1 1 .b L. EVALUATION OF THE REFLECTION ... ....... ... 

2 

Previous results for the range 0S«(> Indicated that thlo section or the problem 
should be recoded using flmpaon’s Rule Instead of the Trapezoidal Rule. It was then de- 
tdod 1 remaining sect tons of the problem aa well, l.e., 

for Uii- parts of the problem in the range -> <«< . Preliminary codes were written and run 
ror those two latter aectlona, yielding satisfactory results. These preliminary codea are 
to re converted to final codes, a relatively minor task, and then run. 


A. Balaer 


ITS B,N. OVERLAP INTEGRALS OF MOLECULAR AND CRYSTAL PHYSICS 
See Section 2.2 of Part I. 


M3 L. 


EIGENVALUE PROBLEM FOR PROPAGATION OF ELECTROMAGNETIC WAVES 


This problem has been do3crlted In detail In Summary Report No. 39. It arose at 
Lincoln Laboratory In connection will the problem of calculating the electromagnetic radio 
frequency field radiated ry a Hertzian dipole to points well beyond the horizon over a 
perfectly reflecting earth through the lower atmosphere idealized to be an Inhomogeneous 
medium with index of refraction decrrns'ng linearly with height. 

Numerical computations were continued during the past quarter. 

H. B. Dwight 
Lincoln Laboratory 

19b B,N AN AUGMENTED PUNE WAVE METHOD AS APPLIED TO SODIUM 

An attempt has been made to weep Input to Its strictest minimum In a program which 
calculates certain electron energies In face-centered and body-centered crystals. 

The program as It now stands makes much uae of punchouts. In the first stage the 
Input consists of the reduced wave vector (specified on Integers) of the point In the 
brlllou n (one at which energies are to be calculated along with the lattice constant and 
the radius or the Inscribed nphere of the unit crystal cell. The output consists of a 
punched tape and an oscilloscope dlnplay of the Information punched. When to this tape la 
added certain other Information, the tape contains all the physical constants neceaaary to 
character ze an augmented plane wave (APW) In ensuing calculations at the particular point 
In the Brlllouln zone. This allows the Input In all other parts of the calculation to be 
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integers. One can regard theae Integers, which are paired, aa the two subscripts which 
characterize an APW. One subscript specifies the order of neighbor In reciprocal space 
to which the APW belongs and the other specifies the "band" to which It belongs. An a 
check, one subscript is always taken to be negative, the other always positive, li is 
then simple to program a check to see If a subscript hsa been missed. 

In the next stage, parameters associated with the APW must be calculated. The 
Input here la the initial tape along with a "subscript tape". The output la a punched 
tape of the calculated parametera followed by a scope display or theae. 

The third stage Is an assembly stage where a certain subnet of APW Is to be 
assembled to form a subspace in which the Hamiltonian la to be diagonalized. r.e first 
part of the Input consists of all the tapeB from the previous stage which contain the compon¬ 
ent APW of the aubapace desired. Theae tapes are read into the compute!- while the program 
checks to see If they refer to the APW at the point required. Thlo feature guards against 
tape mlxups. Then a "oubacrlpt” tape la fed In and completes the input. Aa another safe¬ 
guard, the program checks to aee If all the component APW iaked for nave been read In. In 
this "subscript" tape the order of the subscripts are Important because the secular equation 
routine Tlrat solves the secular equation with one APW. then two, then three, etc., until 
the aubapace la exhausted. The output of this stage la a punched lope which now charact¬ 
erizes a subspace of APW at a point In reciprocal apace. Who! remains to be specified la 
the particular symmetry the APW is to have. 

This latter Information Is the Input of the last stage: generation of matrix 
elements and dlagonallzatlon of the secular equation. The symmetry la specified by rive 
integers. The first of these apeuirtea the dimension or the representation to which the 
APW ure to belong. The olgn of this Integer specifies the behavior under Inversion. 

Three Integers specify the basis of the representation. The I’lrat. two of these Integers 
specify the "projection” operator to be applied to the APW while a fourth Integer specifies 
if the representation or Its adjoint la to be given. If only the dimensionality la given 
the program itself computes the "projection” operator wbloh leads to the fewest number of 
APW In the linear combination of APW which constitutes the symmetriclzed APW. Along with 
these five Integers the tape from the prcvlouo stage completes the Input. Aa output the 
secular equation (energy and overlap matrices) 1 b punched out and displayed, and, of course, 
the eigenvalues and eigenvectors are displayed. 

Much of the routine outlined above has been tested. In the meantime, calculations 
have been carried out at the points (1,0,1 )ir/a, (0,0',0)r/o, and (2,0,Q)i>/o Tor one "band". 
These results Indicate that for all but the symmetric states convergence oven In one "band" 

Is quite good. The convergence for the symmetric states, however, leads one to believe that 
other "bonds" should be Included. Further Information can be found In the quarterly 
Progress Reports of the Solid State and Molecular Theory Group. 

. M. M. Saffren 

Solid State and 
Molecular Theory Group 


195 C. INTESTINAL MOTILITY 

Problem 195 '» 0 study of the effect of radLallon upon the motility or the small 
Intestine In the rabbit. The analysis of records Is being performed using, autocorrelation 
and Fourier transformation, both performed on WVI mder the supervision or D. Hamilton of 
the MIT Servomechanisms Laboratory. 

Analysis of a total or 89 motility records has been completed on WWI Including 
photographic plotting of both autocorrelation and Fourier transform. Seven records ure not, 
yet completed. 

Evaluation of analyzed records has Just Been started. 

Dr. John Farrar 
Evans Memorial 
Hospital 

199 N. LAMINAR BOUNDARY LAYER OF A STEADY, id: PLOW IN THE ENTRANCE REGION OF 

A TUBE 

In connection with the research on heat-transfer coefficients for auperaonto flow 
of air In a round tube, a theoretical Investigation of the characterlnttca of the laminar 
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boundary layar In the entrance region has been carried out. The bourulary layer equatlona 
or continuity, momentum, and energy are to be aolved for apeciric entrance Mach nuabera 
and thermal condltlona at the tube wall. 


Ull'l method la uaed In the numerloal Integration of theae equatlona. 

A new program ualng CS Tor finding the correct Initial condltlona waa run aucceaa- 
fully for the aolutton of the flrat aet or differential equatlona. 

Another program la being prepored Tor the solution or the third aet of differential 

equations. 


203 D, N. 


T.Y.Toong 

Nechanleal Engineering 

RESPONSE OP A FIVE STORY PRAME BUILDING UNDER DYNAMIC LOADING 
See Section 2.2 of Part I. 


20* N. EXCHANGE INTEGRALS BETWEEN REAL SLATER ORBITALS 


Testing on this problem la being continued by P. Merryman of the Physics Depart¬ 
ment of the University of Chicago. 


R.K.Nesbet 
Solid State and 
Molecular Theory Oroup 


212 B,N. DISPERSION CURVES FOR SEISMIC WAVES: MULTILAYERED MEDIA 
See Section 2.2 of Part I. 


215 D. INDUSTRIAL PROCESS CONTROL STUDIES 

Theae studies are being carried out by J.B.Reawlclt and T.P.Goodman of the MIT 
Mechanical Engineering Department under problem numbers 213 and 281. 

The Machine Design Division of the Mechanical Engineering Department Is engaged 
In a project to determine the characteristics of an Industrial process from the random 
fluctuations In the Input end output records during the normal course of operation of the 
process. The method used Is to compute two statistical functions — the autocorrelation 
of the input and the cro8scorrelatlon between the Input and output of the process. Prom 
these correlation functions, the characteristics of the process can be determined much 
more readily than from the original records (ref. 1); henoe the Investigation. These 
computations Involve a large number of successive multiplications and additions, and hence 
can be carried out moat readily on an automatic computer. In thla lnveatlgatlon, correl- 
atlona were performed on Whirlwind I, ualng programs developed at the MIT Servomechanisms 
laboratory. 


The flrat part of this Investigation waa reported In the Quarterly Report Tor 
the last ouerter of 195“ a.*l In reference 2 . 

The Investigation haa been continued In two dlrectlona: 

(1) Operating records were obtained for a dlatlllatlon column and a heat exchanger at the 
oil refinery of the Rock Island Refining Corporation, and correlation functions for theae 
records were computed on WWI. Some preliminary results have been ootalned on the character- 
latlcs or these processes (ref. 5), and some new electronic equipment IS being designed in 
order to obtain additional Information on process characteristics from the correlation 
funct 1 ona. 

(2) As part of a program to extend theae methoda to procesaea with multiple Inputs (ref 3), 
a two-Input process was simulated on an analogue computer and aub.ected to artificial 
random Inputs. Correlation functions for the simulated procees were computed on WWI, 

It la planned to continue thla lnveatlgatlon aa additional operating records from 
Industrial processes become available for analysis. M ^ recoros 
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217 N. VARIATION-PERTURBATION OP ATOMIC WAVE FUNCTION AND ENERGIES 
See Section 2.2 of Part I. 


J. B. Keswick 
T. P. Goodman 
Nochsnlcal Engineering 


218 N. TRANSFORMATION OP INTEGRALS FOR DIATOMIC MOLECULES 

A program has been constructed which converts output Information of Problem 
23“ Into the Torm of Input data required by the programs of Problem IN*. 

This program deals with sets of biquadratic linear runctlonals £ 1 J|lefJ 
which are symmetric under Interchange or Indices 1 and j, of k and l, or or the pairs 
(1J), (kE). These functionals ore actually Integrals with certain additional symmetry 
properties which cause large blocks or them to vanish. The Initial data are the blocks of 
Independent non-vanishing integrals arising from a given set of oasis function:). Prom 
these It is required to construct a tensor of the form 


' l I ^ C f !k I J 4'* I J"] 

The Indices are divided Into a number of ranges, corresponding to basis 
functions of different symmetry types. Indices 1 and J must lie In the same range, which 
Tight be denoted by of, and k and I In a single range A, not necessarily the same as «. 

When ranges*and A are the same they specify a single clock of Integrals £w* | **J ; 
otherwise, there are two, blocks £•■> | * a J and £«#/•< A J from which Integrals must be 

selected In computing AK. In order to locate a given Integral £lj j kfj^on the auxiliary 
drum the Indices must b* J permuted td a canonical order: lij, k * 1 , and (lj)?(ki), 
treated 8S two-digit numbers. The Integrals of a given block occur In this order, which 
ensures that all Independent Integrals are Included and that r.o eciual Integrals are re¬ 
peated. If the index cr has more than one value. Integrals fljlkil'for the various values 
ofcr follow each other on the auxiliary drum. The first oddress of each Integral blook la 
stored In a table In the core memory. 


In the general cape the coefficients C x depend on both the values of « 

and A and their order, so Af . Is symmetric -rider Interchange of liklluea 1 and J ana of k 

and l, but not of the pairs (1J), [id). Por given**, 1, and J 11 5-J) ell elements A^, 

(k t-t) are to be planted In order on the auxiliary drum. The Index A does not run " 
through all ranges of Indices and some may be repeated. However, Index A rune through the 
same ranges In order for all values of «, and** runs through all ranges without repetition. 

In calculations for which thla program is to be used, single blocks [lj I kij "" 
can be too long to be Included In the core memory. Since elements must oe picked out. from 
arbitrary places In the block It appears to be Impossible to avoid making drum tranarerb 
of single CS numbers (or sets of two or three when cr has several values) Instead of the 
relatively more efficient block transfers. The tensor elements A" are transferred to the 
drum in blocks whenever s certain number have been oelculated. Kedn though the program Is 
rorced to use single drum transfers, the over-all speed seems atlll to be limited primarily 
by the CS arithmetic and by the complicated WWI program needed to control the calculation 
ar.d to compute drum addresses. 


R. K. Noaoet 
Solid State and 
Molecular Theory 
Group 
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22A N. COMPUTATIONS OF THE FIELDS OP VERTICAL VELOCITY AND HORIZONTAL DIVERGENCE 

Thin problem determines the field of vertical sir motion Tor a series or weather 
situations. A detailed description of the problem was given In Summary Report No. »l. 

The results will be nnolyzel By the Pressure Change Project under the supervision or 
Proressor James M. Austin of the M.I.T. Meteorology Department. 

The program is running successfully and results are being complied. An unusual 
example oC dnla processing, the program computes and displays about 15,000 flve-dlglt 
words of useful Information on 1A5 frsmeB of film In about four and one-hair minutes. 

W. woir 
Meteorology 


225 B.N. NEUTRON-DEUTEHON SCATTERING 
Sec Section 2.2 of Part I. 


220 D. INVESTIGATION OF THE VORTICITY FIELD IN THE GENERAL CIRCULATION OF THE ATMOSPHERE 
Description of Problem 

The physical processes which are Important In the maintenance or the general cir¬ 
culation or the atmosphere are being Investigated with the aid of a two parameter, non¬ 
linear, quaal-geoatrop.ilc model of atmospheric flow. Tne model makes use of the vortlclty 
equation and the first law of thermodynamics and Incorporates effects of nonadlabatlc 
heating and friction as well as the vertical advectlor. of vortlclty and the transformation 
of horizontal vortlclty Into vertical vortlclty. A spherical coordinate system with 
pressure as the vertical coordinate has been used In the derivation, and the model la being 
applied to the Northern Hemisphere from 15°N to 75°N latitude. The two parameters used to 
describe the now are the 700 mb and 300 mb contour heights. 

With the aid of the model described above It Is planned to examine the role of 
various heating distribution In generating the gross features of the general circulation. 
Ir. particular, solutions for adiabatic flow will be compared with those for: 

a. a simple heating distribution derived from data on the mean monthly heat 
transports across latitude circles; 

b. a heating distribution which takes Into account the so-called "secondary" 
heat aources and sinks which arise because of differences Between land and sea; 

c. a succession or heating distributions of increasing complexity leading to a 
time dependent heating function which varies In the vertical as well as In the horizontal 
direction. 

It la planned to Iterate the solution for a period of several days starting with observed 
Initial conditions. In addition, solutions will be obtained for the case in which the 
Initial conditions are given by a fictitious now pattern characterized by a random dla- 
turbonce 3uperImposed on a basic zonal current. Ti.ese solutions will be Iterated Tor an 
extended period and the day-to-day changes or the flow pattern studied lr. an attempt to 
gain some Insight Into the mechanism by which differential heating generates and maintains 
the observed features of the general circulation. 

The equations of the model are of the following Tonr.. 

(i) °<y> "t ! - N(h,h*, . T( a ( y) 
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* n - -lconst)£’|£ , + J B (h,h>) ♦ const 
JT *• P(h.n-,a m .s T ,P B i *,y) . o 
C) 0(y)!P Z n . H(y) 44 . a^ 

(5) - T - 0(y)* z h- - H(y)^ 

?:,A* a: ‘'‘ ■- ar,: u,e “ill H e coordinates. respectively, F(), 

HI), I(). •*, and .() ar known functions of tne initial conditions, and, 

h - geopotential height of o pressure surface 

h 1 - geopotentlal thickness between two pressure surfaces 

s m* vortlclty at the near, "reasure level 

s^ thermal vortlclty 

P m “ vertical vortlclty 

j.(*,0) • f (y) (i! 22 - 22 ZL ) 

s * l '»a >y a> ay ; 

The following is the procedure Tor solving the equations listed aoove: 

1. Write the equations In finite difference form. 

2. Initially, we are given the fields of h, h*, s , s T . 


3. Solve equation (1) for^ 


by relaxation. 


u . Introduce Into equation (2) and evaluate P ff . 

• 3B 

5. Introduce Into equation (3) and evaluate 


6. Extrapolate jt— and over a short interval of time and add this value 

to the Initial values of h‘ and o to obtain the values of h* and s at the new 

® m 

time. 

7. Use tnese values of h* and a^ In equations (A) and (5) to obtain values of h 
(by relaxation) and s^ at the new time. 

8. Having obtained the new fields or h, h 1 , a ni and s T , Iterate the above pY-o- 
cedure several times until tne desired forecast Is obtained. 

Numerical Procedures 

Wo had planned to solve the finite difference approximations to the oquatlona of 
the model by relaxation, using a graded net of points over the hemisphere. However, we are 
now testing the feasibility of solving by an Iterative application of the lnverae matrlcea 
of the finite dirrerence operatora. Thin method of nolutton appeara worthwhile In view of 
the relatively high apeed of multiplication which Whirlwind I la capable or (approximately 
<•5 mlcroaeconda). Furthermore, this method simplifies the solution for a graded net. 
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Progress ami ('Hart Plana 


Programs ror interpolating the original Cats to obtain values at the grid points 
or the graded net, and for computing the mean and thermal vortlclty have been completed and 
successfully tested. Programs Tor computing -he non-llnesr functions of » and y nave 
teer. written . it no: . *•- tested. A progrer for lnvertlrg a large order matrix (JOth or*0th 
order) by partitioning and using library programs Tor matrix Inversion and multiplication 
(1 th to 20th order) haa Seen written and la currently relng teated. 

Plana for the next quarter Include completion or the matrix Inversions, if this 
me-iod la round to be feasible, and solution cf me equations for at least one initial state. 

Personnel 

Tills research Is being carrU-1 out under the sponsorship o. the Air Porce Cam- 
bridge Research Center by the ruIloM'.ng peraor.a: R.L.Proffer, General Circulation Project, 
’UTi Duane Coale., Oeophyalos Research Directorate; P. Castillo, ueneral Circulation 
Projoct, MtTv Messrs. Cooley ard Castillo are engaged In the programming. 

R.L.Pfeffer 

General Circulation Project 

228 N. EVALUATION OK DIFFERENCE DIFFUSION EQUATION 
Sec Section 2.2 of Part I. 

210 C. DYNAMIC ANALYSIS OP BRIDOES 

Saul Namyet or the MIT Department of Civil end Sanitary Engineering has obtained 
the final results from a long aeries of Whirlwind programs designed to determine the 
dynamic response of a slnglc-degree-or-'reedom system having 8 variety of reslstanee- 
d 1 splice :c- nt functloiin ard subjected to • series of different loading functions. 

The basic differential equation: 
p n - R n • 

la p value ted ly rreana or the aecond order difference equation: 


(dt) x(t ri ) - - 2x(t n ) 4 x(t I1 . 1 ) 


The results of this problem are contained In a project report on Contract No. 
AP 33(bli>)-2208 submitted to Wright Air Development Center, Wrlght-Pstteraon Air Porce 
!«ae, Ohio. 


S. Namyet 

Civil and Sanitary 
Engineering 


2y E.H EIGENVALUES FOR A SPHEROIDAL SQUARE WELL 
See Sect tor ?.i "f Part I. 


236 C. TRANSIENT RESPONSE OF AIRCRAFT STRUCTURES TO AERODYNAMIC HEATINO 

The transient response of alrcrart structures to aerodynamic heating was Initiated 
or. January 3, 1955 by L.A.Schmlt of the MW Aeroelaat lc and Structures Research Laboratory. 
During the second quarter of l'> ‘ . worn has continued on this problem. H. Parecnanlan has 
been responsible for the programming during this report period. 

The over-all problem lo that or Investigating the Influence or aerodynamic heat¬ 
ing on the structural design of high-speed aircraft. One atep In the solution of the over¬ 
all problem la to determine the transient temperature d lstrlbutlona in built-up aircraft 
structure. This phase requires the solution In generalized form of certain Idealized heat 
flow problems. The second Idealized heat Mow problem (Problem III has been formulated and 
solved. (See Summary Report No. 41 for a discussion or Problem I.) 

Problem II can be briefly described aa the thin plate and w«b problem. Heat 
enters the system oy forced convection on the surface of the skin exposed to the external 
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SCHEMATIC REPRESENTATION OF PROBLEM n 


problem la formulated Ir, terms of the following symbols: 

• adlbatlc wall temperature 

- heat transfer coefficient (boundary layer to skin) 

- length of akin considered 

- length of web 

• skin thickness 

• web thickness 

• location along akin 

- location along web 

- time 

- temperature 

■ weight density 
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when I 4 k $ i-l 

(12) K.d+l.K) - | J» w (l.M-a) - [l- f] *»,(!.>«) + * 

when k - ro 


238 B.H. SELP-CONSISTENT CALCULATIONS OP NUCLEAR kail - 

See Section 2.2 of Pert I. 

239 c. GUIDANCE AND CONTROL 


(13) X w UTl.») * jf f» w (l.--l> ♦ [l - jj *,<*•■> 

Once Lx and Ly are assigned, the selection of at In order to prevent divergent 
oscillation or the difference solution Is a routine ratter. (See Reference 1.1 For the 
finite difference equations to be stable the Increments selected Tor Lx, Ly, and Lt must 
be such that the following criteria are aatlafledi 


(16) 

(15) 


1MJ1 i 

Lt 


im! 

LT 
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Satisfactory damping of convergent oscillations has been obtained by trial. 

The rinal production program Is based on a thirty element physical grid (n-m-15). 
Given the numerical values of the four parameters A s , P. 0 and » as well as a value of 
AT the program computes j 9, { from . At the end of each time cycle the complete 

nondlmenalonal temperature distribution |9,} for nondlmenslonal time T ■ lfi” la avail¬ 
able In hlgh-apeed storage. Whenever (1-?,' at J - 0 or (l-9„) at k ■ 0 reaches 0.1, 

O, 2,... ,etc., up to 0.9 at Intervals of 0.1 the nendla.enslonsl temperature dlatrloutlon ft. 
at that time la stored on the magnetic drum. At the end of the computing phase for one 
cose there are eighteen nondlmenslonal temperature distributions, the nondlmenslonal times 
associated with each distribution, and trie numerical values at the given parameters A , 

P, G, X and At. atored on the magnetic drum. It Is possible, with the present program, to 

compute and store the results for as many as ten cases before Introducing the output display 

tapes. The efficiency of the computing phase or this program la good because the basic 

temperature distribution cycle, which constitutes a major portion of the computation. Is 
carried out lr. WVI rather than CS II. 

The first output tape yields a total of seven scope frames per case. The first 
Tour frames plot (1-9) versus x and y for each of eighteen values of t. The fifth frame 
la a plot or (l-9„) versus T for J . 0, 1, 2 etc., up to 13. The sixth frame Is a plot or 
(1-9 ) versus t Tor k « 0, 1, 2, etc., up to 13. The seventh frame Is an alphanumerical 
scope display of the giver', quantities which derine a cane ( 9 a , P, G, .X , and if), the 
eighteen pertinent times for wnlcn temperature distributions are available, and the rinal 
values of f 9jJ • The final values or | would serve as Initial condition In the event 
that lt ever became necessary to carry the calculation nearer to the steady state condition. 

The second output (fp) tape displays all the Information stored on the magnetic 
drum. This Information Is displayed on the scope In decimal fraction form. Since all the 
values of 9 ore between sero and unity the display gives all the values of 9 to four dec¬ 
imal places in numerical, and Immediately usable, form. It should be noted that for Some 
engineering applications the graphical displays provided by the first output tape would 
be sufficiently accurate. 

Work on Problem II la essentially complete. During the next report period,work 
will continue on the over-all problem. The computation of the transient thermal stresses 
which result from the temperature gradients will be considered. 

Reference 

1. Hildebrand, P.B., Methods of Applied Mathematics, Prentice Hall, Inc., New York, 1952, 

P. 323-365. 

L. A.Schalt 
Aeroelaatlc and 
Structures P.csearch 
Laboratory 
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„ Problem was described In tne previous quarterly report \ jmmary Report 

te " 1 nln ?’ th * C:: !! f ,roer ‘" ] ran successfully and Tour runs were 
V °* ore * coGed program (that la, a routine coded In the tlgerralc system des¬ 
cribed in previous reports under Problem 108) has now tee: written for tils problem but 
not tested. 


J.It.turning, Jr. 
Instrumentation 
laboratory 

261 B.N. TRANSIENTS IN DISTILLATION COLUMNS 

The study or an Ideal model of a fractionation tower wiilcn Is undergoing changes 
In lta operating conditions with time has been continued. A program naa beer, written 
which can be used to calculate the plate compositions or a column In wnlch the reed com¬ 
position has undergone a atep change out the product compositions are maintained constant . 
by controlling the llquld-vapor ratios In the upper end lower parts of the column. Pre¬ 
dicted values of llquld-vapor ratios rrorn Instantaneous compositions on plates between 
the feed piste and the top plate are also calculated a a though the column were operating 
at equilibrium. Similarly predicted values of llquld-vapor ratios In the lower part of 
the column are calculated from Instantaneous compositions of pistes between the feed plate 
and the bottom plate. 

Prom these data It la hoped to obtain correlations mat will Indicate: (l) tne 
simplest method for controlling a column by Introducing a delay factor when using a single 
intermediary plate composition In controlling the reflux ratio in the upper part of the 
column or the heat load to the reboller; (2) more exact control of a column by using (a) 
more than one plate composition In each pert of the column, (o) tne time derivative or the 
composition together with Its Instantaneous value, or (e) the instantaneous plate compos¬ 
ition together with the composition of the plate at ace earlier time. 


S.H.Davla, Jr. 
Chemical Engineering 


262 N. NUMBER OF STRUCTURES OF RELATIONS ON FINITE SET 
See Section 2.2 of Part I.- 


265 N. THEORY OP NEUTRON REACTIONS 

In the study of neutron reactions described lr. Su-i.ary Report No. 61, two programs, 
one for finding the power aeries oolutlon to start the solution or the differential equation 
and the other for the solution ltaelf and the cross sections have been written and tested 
separately. They have been combined Into a master program which for sny given x will find 
the cross sections for any range of X‘ by the use of a suitable parameter tape. The program 

has been teated satisfactorily for known values of the cross sections. Trouble of unknown 
origin arose In the first major run but aa soon as this la cleared up lt is planned to run 
the program for / - 2, 5 - .15 flndlng^thc cross sections Tor x - .6(.4) 6.6 and Xg - 

0(5) 160 Tor each xg . stopping when X - . 

The problem Is being programmed for Professor H. Peshbaoh of the Physics Depart¬ 
ment by E. Campbell and E. Mack of the -olnt Computing Croup. 

E. Campbell 

Joint Computing Group 

266 B.N. SCATTERING FROM OXYGEN 

See Section 2.2 of Part I. 
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2U7 C. SURFACE PRESSURE PREDICTION 

This pro: l«r has been carried out under tne supervision or Professor J.G.Bryan or 
the FIT Mathematics Department. The problem, os described In Summary Report No. 41, 
evaluates the Introduction of the quadratic lomponent Into the general linear scheme 
develope 1 • Wadswomti and "alone Tor r .r.-rlcal forecasting. 

The results will he Included In a report to fe submitted to the Geophysical 
Heaearcn Directorate of the U. -. Air Force. 


2 0 C. TRANSLATION PROGRAM FOR THE HUMERICAbLY CONTROLLED MILLING MACHINE 

A preliminary verolon or the NOW translation program naa been completed. 

Teats of thts version have shown that the transition from otralght-llne to circular cuts 
la. If: certain casea. Improperly -omputed ;y the routine. All other featurea described 
in Summary Report No. 41 ore ohtaVnahle with the existing routine, and after this mistake 
la corrected, the syatem will be available In Its entirety. 

The vocabulary accepted By the tranalatlon program has been extended to Include 
a "zero rccd-rate" symbol. A aero feed-rote associated with a cutting Instruction la 
denned to mean that the instruction la actually not to be executed, But that the cut 
preceding or following this Instruction (as appropriate) la to Be modified as If the In¬ 
struction had been obeyed. Thus, the programmer ma;, specify tnat a sequence or cuts Is to 
begin aa If the cutter raid c-)|ce from some point on the work or that It end as If the cutter 
were to go next to some other point. This notation la used to simplify the starting and 
the termination or a continuous sequence or cuts. 

Tir- existing preliminary version has been successfully used by Servomechanisms 
Laboratory personnel Tor programming a few pieces to be cut on the 'CMP. Although the 
oyatem has not yet beer, used sufficiently to permit a reliable statement of lta 
effectiveness, some indication of Its value has .een Obtained. The most difficult of 
the pieces that were tried requires about eight hours of hand computation. Programming 
the same piece for the translation routine required thirty-five minutes, to which must 
be added forty-seven seconds or Whirlwind I computer time for processing the problem. 

The other parts which were programmed Indicate a similar reduction In programming time. 

All of the programs submitted for processing were error-free, despite the unfamiliar new 
notation. 


suit 

more 


The existing error !r t\prellrlnnry version will be corrected, and the re- 
lng routine will re subjected to test. The possibility of extending the system to 
complicated cases, especially to three-dlmenalonal problems, will be considered. 


2h2 N. ANALYSIS OF TWO STORY STEEL FRAME BUILDING 


A. Siegel 
Digital Computer 
Laboratory 


See section 2.2 of Part I. 


?b>' c. A Kill ERA 110} INPUT TRANSLATION PROGRAM 

This two-paas Input translation program was described In some detail In 
Summary Report No. 41. The program translates mnomonlcally coded programs Tor the ERA 
1107 computer punched In Flexowrlter form on paper tape Into the standard bl-octal form 
directly acceptable by an 110? computer. This program will be used by programmers work- 
lng on Problem 120. 

At the moment the program appears to he operating satisfactorily In all re¬ 
spects. However, automatic operation It. the same -.inner as CS and other WWI codlnK 
systems will not be available until the new WWI Input utility program (see Problem 100) 
18 completed. 


J. N. Frankovlch 
Digital Computer 
laboratory 
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E‘7 C. HORIZONTAL STABILIZER ANALYSIS 

Tne problem [a to evaluate me deformations of tne horizontal stall liter of j 
typical present day fighter aircraft subjected to dynamic Input conditions. Although ■ he 
solution or tne problem, has been simplified ny uncoupling the equations for rigid body- 
notions of the aircraft and the equations Tor vibratory motions of tne ata Ulcer, -hi 
rigid oody response solutions serve as part of the Input conditions -o the vibratory 
equations of motion. Two separate programs have Been Bel tp, 
rigid body equations of motion and one for the vibratory equations of mot. Ion. 

The rigid body motions considered ore vertical translation and pitching rols- 
ar -d, therefore, there .re two rigid Body equations. These equations arc lepend.nt 
nonlinear second-order differential equations; the nonllhearlt.- arising In i he ,.-vo- 
dynaalc force terms Introduced ... I . Lnput 'auditions. Because of the non¬ 

linearity , conventional means annot ue used to solve these equations and therefore re- 
-ourae to an approximate numerical scheme la necessary. Tne approximate numerlea] 

-ethod employed !z a numerical Integration representation of each of the unknown varlaLlcs 
vd its derive- '••as In term's of prevloisl; evaluated valuer, or the variable and Its de¬ 
rivative*. Tne acouracj or a scheme of thla nature (step-by-step numerical solution) is 
Impendent on tne increment of time employed In the analysis. The results obtained from 
vnls program for one Input condition nave been cheered satisfactorily at various points 
f. a desk computer solution or tne sure equations. 

To solve for the vibratory motions of the horizontal stabilizer to dynamic In¬ 
put conditions, , dynamic model . presentation of -.he osa and stiffness distributions 
of the Stabilizer wos made. Thla representation consisted of seven lumped masses mon- 


iton of the external forcer. Tne properties oj the spring which represents the stiffness 
llatrlbutlon of the uuekled are.) are altered from their initial pre-budcllng values. 

Aa the oucklc3 continue to lncr"as'-, tne properties of the spring continue to change, 
thereby making the equations of motion nonlinear In form and non-solvable by convention¬ 
al closed form, methods. 

Tne nonlinear equations of the pos:-buc. ling region were linearised and then 
programmed for a step-':.-step numerical solution on th" Whirlwind Computer. The equa¬ 
tions were linearized by dividing the post-buckling region Into parts and making 
straight-line approximations In each of these parts. Although the Individual seta of 
linearized equations "or each part of the poat-buckllns region con be solved In closed 
form, by Laplace Transformations, it Is considerably more expedient to employ a alep- 
1;.-Step numerical solution on the Whirlwind Computer. The atep-by-step numerical solution 
Involves a numerical Integration representation of eacli of the seven variables and Its de¬ 
rivatives In ten s of previously evaluated values of the variable and Its derivatives. 

The program operated satisfactorily for tne seven lumped mass system described 
and for a particular set of Input conditions, it is Intended In the future to analyze 
a twelve lumped masa system subjected to various Input conditions. 

E. S. CrlJelonc 

Aeroclaatlc and 
Structures Research 
Laboratory 

258 C. DYNAMIC ANALYSIS OK AN AIRCRAFT INTEBCEf 

This problem, described In Summary lieoort No. 41, has been completely cheeked 
out, and eight solutions have been run. The first throe solutions were for the purpose 
or determining how well certain force and u-nent equations described the motion or an 
aircraft. These solutions were also used as check solutions for analogue computations 
performed on tne M.I.T. flight simulator or the Dynamic Analysis ana Control laboratory. 

Aa a result of the rirat three solutions and the analogue study, certain 
changes were Indicated as desirable In the given equations and In none constants and In¬ 
put functions involved. These changes were applied to the Whirlwind programs, and new 
solutions were run to study the errect of the changes. 

The results of the second aeries or runs are now being considered by engineers 
of the D.A.C. Laboratory and the Deportment of Aeronautical Engineering. Thraa 
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supplementary solutions are planned to complete the study. 

K. Kavanagh 
Dynamic Analysis 
and Control 
laboratory 

219 L. THE MEDIUM FREQUENCY IONOSPHERIC PROPAGATION STUDY 

The earth's atmosphere between the altitudes o!' approximately 75 and 500 mllea 
Is partially Ionized, and constitutes the region Known as the Ionosphere. Under certain 
circumstances and over certain frequency ranges, the ionosphere Is capable of reflecting 
radio woven; thin Tact la the baaln of most present-day long range communication and 
navigation nystemu. Also, moat studies of.the physical structure of the Ionosphere uti¬ 
lize radio wave reflection as the principal experimental tool. For these reasons, the 
behavior of a radio wave that has undergone Ionospheric reflection la of considerable 
Interest. 


The physical structure or the ionosphere appears to be qulle complicated; It 
Is neither homogeneous, nor Isotropic, nor constant in time. It Is subject to various 
periodic and aperiodic disturbances such as meteors, sunspots, and vagaries in the earth's 
magnetic field. No detailed physical theory of the Ionosphere la Known that satisfactori¬ 
ly accounts for all of the Ionospheric phenomena that have been observed; the problem 
appears to be somewhat worse than giving a detailed theory of the troposphere (lower 
atmosphere) that would satisfactorily Bcnount for the weather. 

However, certain types of Ionospheric problems may prove amenable to a atatla- 
tlcal type of treatment; It la a problem of this type that la the subject of the present 
study. When a steady-s'ate, single-frequency radio wave Is directed at the Ionosphere, 
both the phase and amplitude of the returned wave are observed to vary In time (at a 
fixed point) and in apace (at a fixed time). The statistical structure of this varia¬ 
bility la of Interest. For example, the probability distributions of amplitude arid 
phase, the autocorrelation functions of amplitude and phase, the spatial correlation 
function or amplitude, the cross-correlation function of the amplitude at different 
points of space, and the correlation of amplitude fades and phase shifts are all impor¬ 
tant data. 


In recent years, statistical theories of Ionospheric reflection have been 
developed in which some of the statistical data above are predicted. Some or this theory 
Vs of a simplified nature ( 1 ) ( 2 ) ( 5 ), while some Is Based or. highly refined physical 
theories (a) (s); however, the statistical results are often the same In either case ( 5 ). 
Some of the results predicted by this theory have been applied to Ionospheric studies, 
and have been tested In highly limited ways (b) ( 7 ) (it). However, no thorough study of 
the (short-term) statistics of j reflected wave has yet been made, especially not under 
circumstances appropriate to communication and navigation systems. There exists con¬ 
siderable doubt as to the validity of the theories mentioned. The chief reason why no 
such study has been made in the past Is the formidable amount of computation required; 
the determination of large numbers (several hundred) of correlation functions would be 
Impossible, aa a practical matter. In the absence of large-. y.-speed computers. 


The principal objective of tne present study la to partially close this gap. 
From January through April, 1955, recordings of lonospherlcally reflected waves were made 
lor five hours per day In such a way that the statistical parameters mentioned above 
could be determined. (Under the circumstances of thlo experiment, about one hour la re¬ 
quired lor a statistical sample, so this la roughly equivalent to five independent 
samples per day.) The transmission path was from South Dartmouth, Mass., to Port 
Be 1 vo 1 r, Virginia, a distance (great circle surface distance) of about 385 mllea (620 
Km.). The 1 requency employed was 54j Kilocycles, ana the hours of observation were 
from 1:00 AM to 61 OO AM (local time) dally. Under these circumstances, observed signals 
° r a °' ie - h °P WD'e reflected from what la known aa the E region 


of the Ionosphere. 


The raw data obtained In the experiment hove beer, transcribed onto punched 
paper tape for processing on Whirlwind I, where the statistical parameters ,111 h. ™ 
puted. It Is felt that the results obtained will be of 1nterest P b™h in SaSlc lono! 
sphere theory and in the practical design of communication and navigation sy»?emS? 

_ Thla research was supported from January, lyu, through January Idgft bv the 

Department of State. Since January, 199*, It has been supported Jointly by the’u!sl 


62 


WHIRLWIND CODINO AND AH'LlCATlONS 


Army, Navy, and Air Force. 
Technology. 


Both contracts have been with the Massachusetts Institute of 
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260 N. ELECTRONIC ENERGY OF THE OH MOLECULE 

Most of the work on this problem has so far been concerned with two major parts 
of the calculation, namely, the calculation of the many various integrals Involved, and 
the evaluation of the "Lowdln Overlap Determinants" (LODS). (Ref. l) 

Calculation of Integrals 

The basic set of ore-particle functions are the Hartree-Fock Is, 2 s, and 2p 
atomic oxygen orbitals, and the'hydrogen ground state function. From these, the Slater 
determinants (SDs) having the proper ground state symmetry of the OH molecule (-r) are 
formed. A linear combination of these determinant*], functions is then token as the wave 
function describing the ground states of the molecule, and their coefficients determined 
by the variational principle to give the lowest ground state energy (secular equation). 
The elements of the secular equation are the one- and two-electron interaction Integrals 
between the atomic orbitals. Using a program written for Whirlwind I by P. . CorbatS, 

many of the coulomb, exchange, hybrid, etc., one- and two-center integrals have > een 
calculated so far. Up to several hundred radial Integrations have been computed at one 
time (l.e., in one "run" of Whirlwind). 

The two-center (one- and two-electron) integrals are calculated by first ex¬ 
panding the hydrogen function about the oxygen as center. For this purpose, CorbatS, 
Karo, and Freeman have tested a program for calculating,radial parts 1 9 :' the terms ap¬ 
pearing in the expansion or the functions (ref. 2 ) e-* r and e-kr /r about another 
center. The convergence of the terms in ' 1 'ited thin way seems good. 


Lowdln Overlap Determiner,ts 

Basic to the use of non-orthogoual funatlona are the numerous overlap inte¬ 
grals which appear In the Interaction of SDa. 


A systematic way or handling the extraordinarily large number of tons* that 
arise In this was recently given by Lowdln. An Lowdli. has Shown: the nnn-orlnogons 1 1 ty 
Integral or two Slater determinants !o equal to the determinant of all the non- 
crthSgonallty" integral. ro^.dlYom th,. b.alc sot 0 fcn.-.Uc 
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provide a straightforward way of calculating the matrix elements of the Hamiltonian (or_ 

any other operator) between SDa wade up from a non-orthogonal act of one-particle Tunctlona. 
Even so, the evaluation of the numerous minors and double minors is tedious even for our 
relatively simple case or 9 x 9 SDa. The simplifying features apparent In the use of ortho¬ 
gonal functions are missing, with tne result that many Integrals are not zero in our case. 
Even If ar. Individual LOB is zero. Its double minors are not, and the corresponding matrix 
elements of energy (consisting of two-electron Integrals)glve plentiful contributions to 
each term In the secular equation . 

Lowdln's complete formalism Is being used In our configurational Interaction 
approach to the Oil problem. Thus far the myriad of minora and double-mlnora arlalng from 
the 91 Independent LODa formed from the thirteen baalc SDa have been evaluated. It still 
remains, aa a major task, to have them double checked. 

This work will be used by A.J.Freeman In a PhD thesis In Solid State Physics. 

References 
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261 C. FOURIER SYNTHESIS FOR CRYSTAL STRUCTURES 

During the spring or 1955 several programs were act up for the uae of Whirl¬ 
wind I for crystal structure analysis as carried out In the Crystallographic laboratory 
of the Department of Geology and Geophysics under the supervision of Professor M. J. 
Buerger. 

The principal computational problems arising In crystal structure analysis can 
be described aa follows: 

1. Transformation of the experimentally found "data” (diffracted x-ray am¬ 
plitudes) which are given In "reciprocal apace" to the vector or electron density 
function In "real or crystal" space by means of a Fourier transformation according to 
the general formulas: 


(1) 

f(xyz) -2.2 2 F 


h k 1 

and 


(2) 

p(xyz) - Z 2 1 F 

hkl 


hkl 


cos 2m(hx + ky + lz) 


hkl 


e -2ml(hx + ky + lz ) 


where P(xyz) Is the three dimensional vector density or Pstterson function,/ 0 (xyz) 
the three dimensional electron density function, F . la the diffracted x-ray amplitude 
(in general a complex quantity), p2- F . F* the nltl diffracted x-ray Intensity 
(always a real number), xyz are the coordinates In "real" apace, and hkl the coordi¬ 
nates In reciprocal" space. 

n. .. ?• The second problem Is the reversed process, namely the transformation rrom 
direct to reciprocal space. This amounts to the computation of the scattered am- 

thi formula: “ F hkl * ° f * slmllar Fourl ” transform, according to 


(3) 
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ere r } uelug the scattering power of the J th atom, * y z the coordinates of the 
J ‘ atom and h, k, l the coordinates in "reciprocal" space. ' 


were erene^d v Vhlrll ! lnd iPrograms for special cases of equations (1) and ( 2 ) 

were prepared by 3. N. Simpson, ,r„ Assistant Professor In tne Department or Geology 
and Geophysics and for equation (J) ty 0. Mahoney of tllf WMr> ■ puter staff. 


For both types the moat general (noncentrosymmetrlcal and three dimensional) 
caaea are in preparation. 


These programs were applied by M. J, Buerger to the determine! Ion of the oryn- 
tal structures of Wollastonlte CaSlO, and Pectollte CaN»Jl,0,,H and by Theodor Hahn to the 
computation of the three dimensional Patterson function of dlglyclnhydrochlorlde, 

Trom which It la hoped to obtain the solution of the structure t.y means or the "minimum 
function. 


The Investigations are being continued, and If successful, publication In a 
scientific Journal la envisaged. 


Prof. M. J. Buerger 
Oeology and 
Geophysics 

263 C. FLIGHT PATH OF IN IRC PA FT DURING PULLUP 

This problem was somewhat described In the Summary Report i.o. 111. It In¬ 
volves a set of r simultaneous first order differential equations plus a set of asso¬ 
ciated algeura lc equations. The program for the problem was written for CS I! using 
lhe Horary of Subroutines for sine, cosine, square root, and differential equations. 

The problem concerns high angle pull-up bontlng where the flight path prior to 
pull-up Is a straight-line dive at the target. A pull-up la then executed and calcula¬ 
tions continue until a so-called release point occurs. Tne CS II problem solves the 
oq-atlons Tor each point (At - .1 sec) and considers each point a possible release 
point. A miss is computed Tor each point and oy plotting •Isaea against time the re¬ 
lease points can oe read off. 

This Is done for a variety or Initial conditions and two sets of differential 

equations. 


The programming for this problem was done by C. Block of the M.I.T. Instru¬ 
mentation Laboratory. 


C. Block 

Instrumentation 

Laboratory 


26A C. OPTIMIZATION OF AIRCRAFT ALTERNATOR REGULATING SYSTEM 

In the past, the problem of electrical machine design has been solved by some 
theory, many rules of thumb, and a copious amount of Intuition and Insight. Tne solu¬ 
tion to the proolem by this method may satisfy the design specifications, but the de- 
3lgnen has no reason to expect hlo solution to be an optimum one. By optimum we mean 
In the sense or minimum weight, size, cost, or a combination of size, weight, and cost. 
The amount of calculation required to find such an optimum for a relatively intricate 
machine would be prohibitive If done by hand. 

With the advent or the hlgn-apeed digital computer, the possibility exists of 
rinding an optimum design for an electrical machine quickly and with economy. Tho tech¬ 
niques may be extended to optimize, not only the machine, but a system Including a 

machine. 


Two approaches toward the optimization have been proposed! 

1. Trial and error: An electrical machine may be represented for digital 
computation as a program which determines quantities like heat rate, wave-form distor¬ 
tion, and short-circuit current rrom the dimensions and other physical parameters of 
the machine. The calculated quantities have constraints Imposed on them by the design 
specifications. Using the trial and error approach, we would try many reasonable 
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combination!! or dimensions and parameters, and select tr.e one which satisfies the design 
specifications and gives the optimum machine (minimum weight, for example). This pro¬ 
cedure con be Improved by choosing the changes for the next trial on the basis or the 
results of previous trials - an educated trial and error approach. 

2. As a minimization problem: The optimization problem, mathematically. Is 
one of finding the minima oT a function of - varlacles subject to n constraints. The 
variables are the dimensions ana parameters of the machine, the constraints are the de¬ 
sign specifications, and the runctlon to .e minimized Is weight or size. 

The Tlrat method lias teen used on a small scale using punched card techniques 1 
with good results. home work has been done oy the writer or. planning numerical proce¬ 
dures Tor the second method. These will ce presented In a future report. 

Progreaa In First quarto:- 

The first task that was undertaken waa developing a sufficiently accurate 
calculation scheme Tor use with the computer over the ranges of the design parameters 
to be tested and to compare the reaulta with the behavior of an actual machine. A 
bO-kva aircraft alternator is of particular Interest to this project, and it was chosen 
as the machine to optimize. 

The only special numerical procedure necessary here was s suitable representa¬ 
tion of nonlinear functions, ouch as magnetization curves, for the computation. A third 
order Interpolation routine w/.s written which uses Gausses forward formula.* 

A program has been successfully run to check which design parameters are moat 
Important In matching each of the specifications. 

Work Is In progress to develop a more accurate calculation scheme for the al¬ 
ternator to be used with the optimization procedures. Plano call for the completion of 
the more accurate calculation scheme and tests of the optimization procedure by the end 
of the summer. 

Personnel 

i . sunders - Associate Professor of Electrical Engineering (Vlalting) 

R. F. Nease - Rcaearch Assistant In Electrical Engineering 

J. B. Dennis - Research Assistant In Electrical Engineering 

The work reported here Is being done at the Servomechanisms Laboratory under 
the sponsorship of the U. S. Air Force. 
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265 L. ELECTRON DIFFUSION IN AN ELECTROMAGNETIC FIELD 

Solutions arc desired Tor the second order boundary-value problem 

n m < * eos 2 »y «• gjAy) - g* 2 <y) 


if ly -.6 ' if ly-1.0 * 0 

for values of the parameters £ , , g. 


cm 
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. * r°“ tln * P r0 Ersm was written which used an eat 1 sate 01 *(. I g- t . 

to this Initial value problem. Tnla solution then tt«8 ua< Lf| ettlf tu t 

*(•5) «>» •" ««>unt depending on the degree by w, lab -A | ,, 

solution for g - 2.0 x 10' 5 , 6 . }.6 x 10’ 5 , 1 . • x IB' 5 was determine,:. Changing 

t0 * 10 5 caused the problem to be so sensitive . varlatloi . . , ■ 

variation of .005 In this value caused a varlatloi, ,.f ipproalnaM i . 1 ' • ,, »„ Iv¬ 

or »(y) .6£y £1.0. 


After trying several variations or t-.ls procedure,.. i.jui. if *fe.- ii.ock* 

was abandoned completely In favor of systematic trial - ■ r -. : wclng r.ly tie Information 

thit 3y I y-l.o " as P°* tttve or negative (wit. reaped! • ■ threshold Interval tbout ), 

A subsequent solution has been obtained :y this method. 


D.H.Arden 
Digital Computer 
I-i bora tory 


266 A. CALCULATIONS FOR THE MIT REACTOR 

An Investigation of the response or water-cooled reactor.-, with MTU type fuel 

elements, to slow changes In reactivity la under way. 

A atudy will be made of the resulting time history of reactor power, fuel 
element temperature, and coolant channel pressure. 

The system can be described by three Interacting sets or equations: 1. the 
nuclear kinetics, 2. the thermodynamics, end the fluid mechanics. 

From the thermodynamics and fluid mechanics tne average percentage of steam in 
the coolant channela Is calculated 88 a function or the heat output of the renotor ami the 

coolant flow rate. 


Assuming fog flow, the properties of the steam-water mixture can be calculated 
and used In the Fanning equation. 

(1) 


? n * P m * 


i! f, v _ v 1, AD 

gS 2pmV' 2r, J* 


An energy balance between the entrance to the coolant cnanncl and point i can be written 

hZ„ „ .,2 L. 


( 2 ) 


W T 0> * -v 2 - a - jgr <•* - •*> - ? 


These two equations, together with tabulated date, con be solved by trial and 
error for the pressure drop In the fuel element and the percent of "team. 


(3) 


For the nuclear kinetics the following equation Is assumed: 

1 

1 






fission 


4 Vi * ' g 


The reactivity can be written as 

<*) P-Pe+Pm 

in which P e is the externally controlled reactivity and p m changes with tne amount of 

steam In the reactor. 

The calculation of p m la described below. 

The critical maoa of a critical reactor, or the period > for non-orltlcal 
reactor can be found by solving the equations: 
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<R) D s V - •F-’Ps + pVi'O 

o 

(O + <*• D f* - 

The two-group model was assumed. 

— was assumed - 0. 

r r 

Only the highest eigenvalue solution of A is considered to give 0 - 0 Q e t . 

For the calculation of the critical maos > - 0 and the solution Is obtained 
by varying the critical mass. 

For the calculation of the period, the critical mas3 la constant and the 
solution ts obtained by varying the period a . 

The solution Is obtained by writing out the appropriate differential equations 
fir the fuel and reflector regions and assuming the following conditions: 1.) the fluxes 
are continuous at the boundaries; 2.) the net currents are continuous at the boundaries; 
5.) the flux is 0 at the outer boundary; A.) spherical symmetry Is assumed. 

In this proolen a four region reactor 18 calculated: non-multlplylng region 
In the center, furl In the next region and two non-multlplylng regions on tiie outside. 

In order to have a solvable set of equations, the determinant of the coefficients 
of the resulting 12 equations must equal zero. 

The correct value of \ or the critical mass Is found by trial and error. 

In order to bypass difficulties In evaluating determinants with values close 
to zero, the following procedure was followed. 

One of khe arbitrary constants was assumed equal to one. Eleven of the twelve 
simultaneous equations were solved and the results substituted In the twelfth equation. 

It was postulated that If the 12th equation were close to zero, all the other 
equations would oe close to zero. 

This turned out to be true In all cases calculated. 

with the aid of perturbation theory It was possible to obtain expressions for 
the average thermal neutron lifetime, f (Ref. 1) and the fractional change In k per unit 
fraction of steam produced in the fuel region, also called void coefficient (Re?? l). 

-p 5 "* V? ?s dV 



- p *r V^r‘ lv ' Z [n 2 o), ^s' lv - * f -I \/ a / v g*S ^dv-?*,.^/ V0 , ^ 

/ v 0> 3 dV 


k i. 
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These two quantities give Insight Into the behavior o' /, re,dor In n transient. 
The. will be calculated In earn: case studied. Knowing the pump characterin' ics it will he 
P?s 3 l"lc to find Tor each power level the flow rote and the percent i team. This in turn 
will give the change In reactivity calculated from perturbation theory. 

Finally, after these values have seen calculated and tabulated on wwi, a solution 
of equation ( 5 ) will e nought on an analogue computer. 

Tl.c thormodynamlca unn hydraulic calculations were net up ;-y h. Howard, Dynamic 
Analyul: ana Control Uicorator,. The nui ill programmlra roi 

•1 Mfi doni . Troost. This prog, , carried oul forth' 

Project. Ho definite results have been 0 1 

Homeric la tui'e 

1 ’ • pressure (lc/ft' ) 

W - total flow rote (lu/sec) 

v - specific volume (ft'/lb) 

- length (ft) 

* fraction factor 

* V L n < rt > 

- hydraulic radius (ft) 

- 52.2 ft. lbs. matter/aer. 1 .force 
» area for flow (ft ? ) 

- temperature (°F) 

= heat of vaporization (yr^) 

- steam flow rate (lbs/sec) 

= total heat added In section considered (BTU/sec) 


L 

r 

AL 

r h 

g 

S 

T 


Cl 

J 

C P 

subscripts m and n refer to 2 different points nlong the length of the fuel element 
subscript o refers to the entrance 


constant ( ’ t T lb - ) 

BTU n 

apeclflc heat (DTU/lb F) 


e 

? 

V 

2 , 

°«i 

c, 


t 

7 

D 

0 


» reactivity 

«. Traction of delayed neutrons 
» average number of neutrons per rioalon 

- macroscopic fission croas-oectlon 

„ decay constant of l’th group of delayed neutron precursors 

- concentration of 1 ’th group of delayed neutron procursors 

- mean thermal neutron lifetime 

- time (seconds) 

» the voctor operator del. (or naola) 

- diffusion coefficient 

- neutron flux 
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1 

k 

P 

r 

A 

v 

aubsorlpt 




• macroscopic cross-section 

- Infinite reproduction factor 

- resonance escape probability 
■ average speed of neutrons 

- period of the reactor 

• volume of the reactor 
s la for the slow group 
f la for the fast group 

tr la for transport cross-section 

• Indicates an adjoint flux 




Reference 

1 . The formulas used were developed from perturbation theory by Dr. M. Benedict 

M. Trooat 

Chemical Engineering 

367 B. NUMERICALLY CONTROLLED MILLING MACHINE TURBINE BLADE 

Machining warped surfaces on the NCMM Involves the determination of a set of 
curves, each of which In to be approximated by a series of straight line segments. These 
curves are the Intersections of the surfaces with a series of plane3. The straight line 
segments used to approximate each curve lie In the plane whose Intersection with the surface 
forms that curve. The path of the center of a cutting tool during the machinery operation 
In a series of straight lines between points on the normals to the surface and at a dis¬ 
tance Trom the surface equal to the radius or a ball-shaped cutting tool. 

The given data locates certain points on the surfaces of a turbine bucket for an 
aircraft Jet engine currently In mass production. These points are given for the boundaries 
ol' eight croBs sections on the airfoil portion of the bucket. At present, the problem Is 
one of Interpolation between given points on the cross sections and also between cross 
sections In order to define the surfaces In greater detail. The Interpolation will result 
In a series or points representing curves formed by the Intersections of the surfaces with 
planes parallel to the given cross sections. 


A program to perform this Interpolation nas been written and Is now being checked. 
It is bssed on the formulo for polynomial Interpolation In terms of divided differences 

fix) - y„ ♦ [x 0 Xj] (x-x 0 > -t [x 0 x,x 2 ] (x-XqHx-Xj) * [x^x^J 

(*-* 0 )(x—Xj)(x-Xg) + .+ [w-- x n] (x-x 0 )(x-x 1 )...(x-x n _ 1 ). 

Due to the error Innerent In taking measurements of »he giver, data from experimental models 
of ali'ToIl aecUons, Divided differences of higher than the third order show increasing 
divergence. The above formula was applied, therefore, to the third order divided differ¬ 
ences. Jork In the Immediate ruturc will consist or calculating the normals to the given 
surface at each point found in the above program, locating with reference to a set of 
orthogonal axes points on each normal at a distance from the surface equal to the radius 
?'A tInE e ? ol » programming these points as the path or the tool center for the 
■ICMM, At present. It Is planned that this last program will oe run on the NCMM and that 
the airfoil section of a pair of forging die blocks will actually be machined. 


or be ? U " w ’? h , lh * ltlea of Mltlne a Called economic evaluation 

or the conventional methods of machining the dies In question as compared to the method 
outlined above, ulnce this project Is for academic credit toward a Bachelor's degree, time 
Umltations have forced the elimination of further study Into the economics or the two 
methods. 


G. Bromfleld 
Business and 
Engineering 
Administration 
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?6B C. EXTRAPOLATION TECHNIQUE: 

The proDlem under conalderatlor. waa the extrapolation or certain time rerlei 
arising from a particular, classified physical problem. Physical considerations made 
acceptable the hypothesis that the true series could ne represented v means of low de¬ 
gree polynomials. 

The following two methods of determining those polynomials were considered: 

1 . find the polynomial which Is neat In a weighted, least squares sense over a 
certain set of data points, or 

3. first compute n+l equally spaced and anoethed data points and then pass a 

polynomial of degree n througn these points. 

Whirlwind I was used to apply these methods to the particular data under con- 
alderatlon. From the results of these runs ■ -len of the two methods tiara com¬ 

pared In terms of the parameters of each method. 

The second method proved to be superior for the data used. However, no general 
conclusions can be deduced from the study -vide. 

: . A. Petrlck 

Scrvomeohnnlsms 

Uuorulory 


260 D. DYNAMIC BEHAVIOR OF SHEAR 

See 3ectlon 2.2 of Part I. 


270 B. CRITICAL NASS CM/’-. '.LINDAICAL GEOMETRY 

There are no exact analytic methods for calculating the critical masses of re¬ 
flected reactors or nor" tnar. one variable il-enslon. Spherically symmetric reactors 
can ve calculated exactly, but cylindrical and parallelepiped ones cannot. One Is forced 
cither to a multlgroup-mul-lreglon point-ny-point Integration, or to finding some approxi¬ 
mation to the exact solution. : 1 a problem will it tempt to test approximations for cylin¬ 
drical geometry. 

Experimental critical -asses of cyllndrlcnl reactors were obtained by A. II. 

Snell In 19A6, for core- containing uranium 236, heavy water, and aluminum, surrounded by 
a heavy water reflector. If a'satisfactory approximation were developed. It could be e»- 
tended to the designs of other reactors, notably the projected MI - " rear'or. 

The approximation now being rested Is patterned artor one developed oy J. F. 
Hill. The exact analytic solutions are written for two coses: a radially in I 
reactor and an er.d reflected one. Both solutions are Obtained from the basic neutron 
diffusion equation, a 2 ?, . * "lS*, , - 0. Two groups are used (fast n.-.d thermal) and two 
regions (core and reflector). The boundary conditions arc that fluxes sn- neutron cur¬ 
rents be continuous M tne core-reflector Interface, and that the rluxea vanish at the 
outer surface of the reflector. 

Applying boundary conditions, the solutions are of the form (Tor the radlolly 

reflected case) 


AX * CY - F.'j 


CY'] . B lr FV 
FjAX * EjCY - SjF/j t 07 . ? 

D le[-V X ' * S 3 CY j • U lr ['3 F V * °V] 


Fast Flux 

Frst Current 

Slow Flux 

Slow Current 
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By Cramer'a rule, the determinant of the coerriclents, A, C, F, and G must equal 
aero Tor non-trlvlal aolutlona. Similar equations are found for the end reflected case 
except that X, V, and Z are trigonometric and hyperbolic functions. 

The aolutlona are then superimposed. One obtains (In Snell's case) a fourth 
and eighth order determinant which must simultaneously equal zero. All physical con¬ 
stants and dimensions are fixed, leaving two variables, critical mass and the radlal- 
longltudlnal Duckling split, which correspond to the two determinants to be solved. 

Progress to date has been the writing of the program, and the development of a 
subroutine to calculate the Bessel functions, J Q , Jj, I Q , Ij, K Q , and Kj. No critical 
mass calculation has yet been obtained. 

An additional section has been written, but not tested, for calculation of the 
neutron lifetime and void coefficient. This la accomplished by taking a 100 x 100 mesh 
net covering the reactor, and Integrating the fluxes numerically. First order perturba¬ 
tion theory Is used. 

This technique will also be useful for calculation of reactivity losses due to 
burnout poisoning, etc. In addition, plots of the radial and longitudinal fast and slow 
(real and adjoint) fluxes can be obtained. 

If the method proves feasible. It will be extended to a four-region reactor, 
possibly with more energy groups. 

Another approximation method exists, and it Is planned to try It In the future. 
It Involves expanding the fluxes In Pourler-slne and Pourler-Bessel series, and solving 
the resulting equations. It has the advantage of using three energy groups, and should 

be more accurate. 

Nomenclature 

Dj e - Fast core diffusion coefficient 

Dj r - Fast reflector diffusion coefficient 

D 2c * Slow core diffusion coefficient 

D Zp * Slow reflector diffusion coefficient 

R c • Radius of 'cor* 

R_ • Radius of reflector 
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Sj, S,,, 5j - Coupling coefficients 

X - 3 o 0Vc) 


x ofrr»o) 


T o^ir R c) 


! 0 ftlr B r1 
M*lr n r) 


K o(*lrM 


• Radial part of real buckling solution A* 2 

• Radial port of Imaginary buckling solutionV 2 
■ Analogue of reciprocal fast dirfualon length 
- Analogue of reciprocal slow diffusion length 


*r\ -/•* 




tir-iV 'A 

t!r-{fc v*-/*5 

T * 8 6 e In reflector 

L S - tncrmal diffusion length In reflector 

This study Is being carried out In connection with an Sc. D. thesis for the 
Chemical Engineering Department. 

J.H.Powell, Jr. 

Chemical Engineering 

271 B. EVALUATION OF A BEAN SPLITTING TECHNIQUE 

This problem Is a study of a proposed method of Improving the accuracy with which 
conventional search radar systems determine target azimuth angles with respect to a rixed 
coordinate system. The operation of extracting aslmuth Information with an error leas than 
tne radar antenna's beam width la sometimes referred to as 'beam splitting'. The necessary 
programming details for WWI are being performed by Peter F. Engel, and the numerical results 
Obtained will be Incorporated Into an K.b. thesis report to ue submitted by him to the 
Department of Electrical Engineering at FIT In August, 195'j. The problem represents one 
pt-ase of a comparative analysis being conducted by Mr. Charlton K. Walter of the Applied 
Mathematics Section of the Computer Laboratory of the Air Force Cambridge Research Center. 

Vr. Robert Bernstein of the Electronics Research Laboratories at Columbia University was 
responsible for the application or a maximum likelihood estimation procedure to the beam 
splitting problem, and hla work la described In detail In a report (hcf. 1) published by 
Columbia University. 

The beam splitting problem Is characterised. In part, by Inherent statistical pro¬ 
perties, and the particular solution under consideration here Is essentially an application 
of the maximum likelihood estimation of a population parameter. A principal objective Is to 
compute estimates of a target's azimuth position using simulated radar data, and to study 
the dispersion or these estimates about the known position. 

The type of search radar system considered employs a continuously rotating antenna, 
radiating one or more directive beams of pulsed signals on an S-band microwave carrier. 

Pulse widths are In the vicinity of one tentr to four microseconds and are emitted .it the 
common pulse repetition frequencies of 200 to boo pulses per second. The radiation pattern 
of the antenna has the form j sln Y j? where y is the angular deviation from the antenna axis. 

As the radiated beam Illuminates a target, reflected pulses are returned to the 
receiver over an arc length equal to the antenna's beam width. IT there were no random 
noise or target scintillation present, the pulse amplitudes would follow the antenna pattern. 
Ir. general, however, noise and scintillation combine with the echo signal and produce a 
skewed or distorted pattern from which It Is difficult to determine the target's azimuth 
position by inspection. By considering each pulse In the azimuth scan as a single sample 
of a random variable whose probability distribution Is known, the target azimuth, 0 , can be 
estimated. 

A threshold, h. Is set arbitrarily to quantize the Individual pulse* making up 
s complete scan Into two discrete levels. The probability of a pulse being largo enough to 
exceed the threshold level depends upon the gain of the antenna In the direction or the 
target at the Instant the pulse was received. Thin In turn depends on tne angle between the 
antenna axis and the azimuth of the terget at that lnatont. The azimuth of the antenna's 
axis Is always known, thus It Is possible to formulate the probability, Tj (x,;0,K), of 
receiving a pulse of amplitude Xj, as a function of target azimuth 0, and the algnal-to- 
nolse ratio, K. If the simplifying assumption la made that successive pulses are statis¬ 
tically Independent, a Joint probability, F(xj....x n ;8,K), for detecting all the pulses of 
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a scan can be formed. This distribution la called the likelihood runctlon, L, for that 
particular scan. 


WHIRLWIND CODINO AND APPLICATION,': 


L - F(x,.*„ i 8,K) • * i «.*> 


where n la the number of pulses per acan. The maximum likelihood method now conalata In 
choosing as an estimate of the unknown population value of 9, the particular‘JS*- 
renders L ns great us possible. .'Wee log L attains lta maximum for the same value or 9. 
as L, one haa to solve the equation, 


a ioft^L _ £ 3log fj (x, 1 9,K) 


The solution, 3, Is the desired target azimuth. Equation (2) can be reduced to a more 
convenient Torrr, 

(?) w^u.ajx, - o 

where W (K,0) la a weighting function of the two variables K and 9. For computational 
purposes, V, can be plotted as a runctlon of a. Tor discrete values of algnal-to-nolse 
ratio K. 1 

After determining 8 Tor each scan, the computer compares the answer with the 
known position of the target, and records the difference between the two. This la re¬ 
peated for several video threshold levels, h, and for several values or r., to study the 
distributions of the dispersion In 8. 

Runs have been completed during which position estimates were computed at four 
discrete thresholds. In all four cases, tot azimuth scans were processed, each of which 
contained 35 pulses. The resulting distributions of the estimates were found to have 
their mean values displaced from the true position by not more than 3 o/o of the antenna's 
beam width, and each distribution had a standard deviation about the mean of approxi¬ 
mately 9 o/o or the beam width. The antenna's beam width Is taken to be the distance be¬ 
tween the 3 db points of the radiation pattern, and for the results cited above, was 
equal to one degree.^ 

References 

1. Technical Report 7-1/128; Department of Electrical Engineering; Electronics Research 
laboratories; Columbia University Engineering Center; b32 West 125th Street, New York 27, 

New York. 

P. F. Engel 

Electrical 

Engineering 

272 L. GENERAL RAYDIST SOLUTION 

A system of equations has been derived by Dr. G. C. Sponsler of Lincoln Labora¬ 
tory to provide for a particular aet of data the position of an aircraft In latitude and 
longitude to an accuracy of better than alx parta In one million. 

The Raydlst system employs one master radio receiver and three slave stations, 
any two of which may be used with the master station to provide data for 3 solution. A 
system or rectangular axea la referred to an origin located on a plane determined Dy tne 
master and two slave stations. The aircraft Is assumed to be on an ellipsoid. A system 
of three simultaneous quadratic equations Is then derived, the coefficients of which are 
functions of Raydlet data, and the solutions of which are rectangular coordinates of the 
elroraft. These rectangular coordinates are then transformed with latitude and longitude. 

A program was written to solve the simultaneous system of equations by minimizing 
a sum or squares. It wss learned, however, that convergence to the solution was too slow 
to render this method practical. Therefore, the system of equations was reduced to a single 
eighth degree polynomial equation by eliminating variables. A program for computing the 


coefficients of this polynomial and obtaining lta roote la now being written. 


0. Mahoney 
Digital Computer 
laboratory 


ANALYSIS OF AIR SHOWER DATA 


or wit i. d ; rcc “? n of Professor Bruno Roaal or the Fhyalca Department 

or ..IT 18 engaged In an experimental study or the size spectrum and arrival directions or 
large air "lowers. This group includes 0. Clark, J. Earl, w. Krauahaar. J. Linsley, and 
F. Scherb. 0. Clark and F. Scherb are primarily responsible ror the design of tne VWI 

analysis program. 


The data for this study la obtained with twenty large scintillation counters 
arranged In an array 500 m. in diameter. The density of particles and the arrival time 
of the shower front at each of the twenty counters la recorded whenever an air shower 
strikes the array. From these data the size of the ahower, the location of the core, and 
the direction of arrival are determined by digital computation. 

The data consist of particle densities g, (1 • l, 2, ...,20), relative arrival 
times tj, the sidereal time, and various instrumental constants Buch as tne counter locations 
X,,y,. The direction cosines f, m, n or the shower axis relative to the ground are deter¬ 
mined by finding those values of k, l, and m which minimize the expression 

'X- ^ (k-fXj - my j - ctj) ? 

n la related to l ar.d m Oy the equation 
n 1 - f 2 - m 2 

Tne right ascension and declination of the shower axis are then computed rrom the direction 
cosines and the sidereal time. 

The core location (X,Y) end shower size N are found by determining the values of 
, £ , X, and Y which minimize the expression 

*■ !t (f t - 8l > 2 

i-i 1 ■ 1 


B, =/(X - »,) Z ♦ (Y -y,) 2 - [KX-x.)4 mCY-yjlJ 2 

and w ( la a weighting raetor. The expression for f ( Is a convenient empirical represent¬ 
ation of the lateral density distribution or particles In an air shower. The size or the 
shower la then computed according to the equation 


The rirst part of the computation program la devoted to the reduction of the 
experimental data to a fori sultatU for computation. The next part la a solution or the 
three first degree simultaneous equations for /, and n (which come from the minimization 
requirement) and the computation of celestial coordinates. The lest part of the program Is 
the minimization of the expression * by means of Iterative subroutine, and tne calculation 
of N. 


7 * 
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The program ha a been written and parta of It have been teated. 


Effort will now be devoted to teetlng the program. When the program la aucceaa- 
fully run. a thorough atudy will be made of the accuracy and resolution of the experimental 
and analysis aetup ualng flctltloue ahowera. Actual ehower data will then be analyzed. 


O. Clark 

P. Scherb 
Phyalca Department 


27A N. MULTIPLE SCATTERING OF WAVES FROM A SPATIAL ARRAY OP SPHERICAL SCATTERERS 
See Seotlon 2.2 of Part I. 

275 B. BUCKLINO OF SHALLOW ELASTIC SHELLS 
See Section 2.2 of Part I. 

277 C. HORIZONTAL STABILIZER MODES. SHAPES AND FREQUENCIES 

Thla problem lnvolvea the determination of natural mode ahapea and frequenclea 
of an airplane atructure repreaented by a lumped maaa dynamic model. The aet of almultan- 
ooue linear equatlona which governa the rree vlbratlona or the ayatem la aolved by a pro- 
scan or matrix Iteration to yield the mode ahapea and frequenclea. 

The aet of almultaneoua equatlona may be written In matrix form aa: 

MM -jj*l 

The A'a and w'a are the modal coefflclenta and the frequenclea reapectlvely. The D matrix 
la computed by combining the deflection Influence coefflclenta Tor the dynamic model, the 
conatanta pertaining to the maaa and geometric propertlea of the aircraft atructure, and 
the conatanta ror the rigid body degreea of freedom of pitching and vertical tranalatlon. 

The main tape for the program conalata of a aet-up for computing the D matrix. 
Iterating for each mode ahape and frequency, and a proceaa for aweeplng the D matrix to 
eliminate each aucceaalve mode after convergence. The tape ha a been prepared In a general 
form which llmlte neither the alze or the D matrix nor the number of modea to be found. 
However, the alze of the D matrix and the number of modea to be found are preaently limited 
by the hlgh-apeed storage capacity of the computer. 

• 

A secondary tape, which lo fed Into the computer with the main tape, Includes the 
matrix of Influence coefficients, the constants associated with the two rigid body degrees 
of freedom, and a constant equal to the number of modes to be round. If It la so desired. 

It la a simple matter to eliminate either or both of the rigid body degrees or Treedom and 
obtain cantilever mode ahapea and frequencies. 

Three problems have been aolved during thla period with satisfactory results. 
First a cantilevered beam model, consisting of three lumped masses, was run to determine 
the soundness of the program. Next, a cantilever model with seven meases was run and seven 
mode shapes and frequencies were obtainod. Finally, a 12 x 12 matrix containing both pitch¬ 
ing and vertical tranalatlon degrees of freedom successfully yielded the first five mode 
anapea and frequencies. 

It Is contemplated In the future to run investigations on various other models. 

It la also possible to use thla program lo solve other eigenvalue problems which arise In 
other branches of engineering. 

N.P.Hobbs 
K.R.Wetmore 
Aeroelaatlc and 
Structures Research 
Laboratory 

278 N. ENEROY LEVELS OP DIATOMIC HYDRIDES (L1H) 

In e Quarterly Progress Report of the Solid State and Molecular Theory Group at 
MIT (Ref. 1) we considered brleriy the problem or the electronic energy of the lithium 
hydride molecule. Using la and 2s atomic wave functions for the lithium atom, and a le 


- -—-—avuika mow di* bii ico occurred with proper ground 

symmetry and propoeed to examine the problem of Interaction between then* piston as a 
function of the Internuclear distance. 


The specific atomic functions which we are ualng are the orthonormal le. 2a. and 
2p Hartree-Pock wave functions for lithium [R,f. 2), and, located on the hydrogen atom, e 
normalized la Slater Atomic Orbital which la expanded about the lithium atom for purposes 
of calculation. Thla procedure of expanding an orbital located on one center about another 
center has been outlined uy P.J.CorbetS (Ref. J) and le particularly simple for the analytic 
rorm we have ehoaen. The Hartree-Fock wave functions are, of oouroe, orthogonal; but the 
hydrogen wave runctlon has not been made orthogonal to the lithium wave functions ao that 
the overlap lnteg -ala must be evaluated. 


All of the one-electron and two-electron Integrate have been obtained Tor the 
lithium atom taking Into consideration the la, 2a, and 2 p funotlonsi and all of the add¬ 
itional overlap, one-electron, and two-electron Integrals for the lithium hydride molecule 
near the observed Internuclear distance, using only the le end 2 a lithium functions end 
the hydrogen la runctlon, have also been calculated. With the exception of the kinetic 
energy Integrals the computations have been performed on the Whirlwind digital computer us Ire 
the routines developed oy P.J.CorbetS (Ref. 5). 

An Intermediate calculation hea been completed for the electronic energy of the 
neutral lithium atom aa a partial check of our work. Using the le, 2a, end 2p Hartroc-Foek 
lithium wave functions one may construct the possible determlnantal wave functions for the 
end zp cases; from these, four atatea may fie rormed which poaaese the Z S ground state 
symmetry and rive atatea with ‘T symmetry. Ilie single determinant representing the unexcited 
ground atate gave an electronic energy within one-hair per cent of the experimental value, 
In good agreement with that found by Pock and Petraehen (Ref. A). When the matrix of Inter¬ 
action of the four ‘S states was diagonalized, there was an almost Insignificant decrease of 
the electronic energy for the ground atate. This eeems to Indicate that the configuration 
interaction with the excited atatea rcaulta In a negligible correction. The Z P atate waa 
similarly but slightly lowered when configuration Interaction with the four other excited 

S tates was Included. The calculated value for the wave length of the transition between the 
P and Zs atatea waa also within one-half per cent of that experimentally observed. 

A preliminary valence bond calculation haa been completed for the lithium hydride 
molecule near the observed Internuclear distance. The value of the electronic binding 
energy for this single state waa found to be about one-third of the observed value of 
2.5 + 0.2 electron volts given by a linear Blrge-Sponer extrapolation (Ref. 5). 

Attention at the moment la centered on oomblnlng the basic Integrals which we have 
obtained, to form the Hamiltonian matrix and the overlap matrix. We shall bo able then to 
determine the extent of configuration Interaction In refining the valence bond calculation 
which haa been made. We are also'prepared to extend tho treatment to a number of Inter¬ 
nuclear distances In the near future. 
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280 B. CORRELATION FUNCTION 

Thla problem concern, the application of the theory o' - Stationary Time Serial to 
Meteorology. Specifically, It 1. concerned directly with the methodB and theory of JA"*" 
Prediction as developed by Norbert Wiener which require, extensive uae of harmonic analyala. 

The work repreacnta an attempt to oxtend the ao-called factorization to more than 
two aeries and to set up the problem of determining the prediction operator Tor multiple 
aerlea on a practical computational baala. So far as la known, no one hay aa yet obtained 
the prediction operator for two or more aeries by direct application of Wiener a theory. 


horizontal ConnerMone * average redundancy and It waa ale round that networks wi.lc, i.,.d 
Iona ?t lei betw^n memnii; 1' il", . V' TV®" ln * >••*>>) rather Urn.. vertl .,1 n,n„.el- 
tlon >ilwnnn c!^nO n^ r , ?,. dJO '" K 1,w,la > **>‘l»lted an Increasingly non-linear rela¬ 
tion between redundancy and the number of stepa toward completion or information iraneter. 
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The computational work can be separated aa follows: obtaining correlation 
functional obtaining spectra of theso correlation functions and related runctlonas factor¬ 
ization (which leads directly to the prediction operator). At the present time, the first 
phase la eaoentlally completed, it la expected that the entire problem will be completed 
during the next quarter. 

The programming has been performed by R. E. Huachke, and B. Rankin Is aiding In 
the theoretical considerations of the problem. 
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28^ N. AUOMENTED PUNE WAVE METHOD AS APPLIED TO CHROMIUM CRYSTAL 


Thin research ln being undertaken aa a doctoral theals ln the Department of 
Meteorology, Massachusetts Institute of Technology. 

P. Ha nna 
Meteorology 


The radial Schroedlnge:- equation Integration routine haa been modified to display 
wave functions numerically on the scope. TntB facilitates finding the range or energy Tor 
which the wave function remains essentially unchanged at a given radius, with thla infor¬ 
mation, aa haa been pointed out by Leighton, the Integration can be apooded up conalder- 


281 C. CORRELATIONS AND TRANSFORMS 

See report under problem number 215 above. 

282 B. HELICOPTER BLADE FLAPPING INSTABILITY 

Equations describing the behaviour of a helicopter rotor blade having three 
degrees of freedom, flapping, bending, and toralon, were formulated. The analysis, which 
took full account of reverse flow effects, led to three differential equations each of 
which contained a term Involving an Integration over the span of the blade. These equat¬ 
ions were solved by Gill's method, a step-by-step numerical integration procedure similar 
to the more familiar Runge-Kutta method. A program making use of this method was drawn up 
for use on Whirlwind I, 06 II. The program will determine, for a particular flight con¬ 
dition, the transient response of the blade by evaluating and recording on the oscilloscope 
values of each of the degrees of freedom at fifteen degree Intervals of the azimuth angle. 

It was found that the computer required approximately eight seconds to compute the values 
at each azimuth angle. 

Runs were made to determine the effect on blade stability of a variation in each 
of the following: 1 ) reverse flow, 2 ) blade weight, }) pitch-flap coupling angle, 4 ) rotor 
angular velocity, and 5) Inclusion or a toralon degree of freedom. 

The conclusions resulting from the above were: 1 ) inclusion of reverse flow effects 
lead to higher values or forward speed at which Instability occurs than would otherwise be 
found, 2 ) at the night speed Investigated, the heavier the blade, the more stable It Is, 

5) on the basis of the runs made (only two ln number), the pitch-flap coupling angle does 
not. appreciably affect the damping of the flapping motion, and 4) due to lack of time, the 
author was not able to obtain any results from A) and 5 ) above. However, initial results 
obtained while the bugs' were being workod out of the tapes seem to indicate that the tor¬ 
sion degree of rreedom has an Important destabilizing effect. 

A complete description or this problem, including all Whirlwind results, la pre¬ 
sented In a Master's thesis submitted to the Department of Aeronautical Engineering. 


Since all the routines described under problem I 9 I 1 ore used for thla problem, 
some of those routines have been tested under this problem number. 

M.M.Saffren 
Solid State and 
Molecular Theory 
Group 

286 B. RESPONSE OF THE HUMAN PILOT IN A DAY SUPERIORITY TYPE FIGHTER FLYING A LEAD 

PURSUIT COURSE 

This problem was undertaken for partial Master's thesis credit by Donald 
Spangenberg and William Banka. The programming was cone by Charles Block of the MIT 
Instrumentation laboratory using the "George" code (.-in algebraic system described ln 
previous Summary Reports under Problem 108). It Involved solving 25 simultaneous first- 
order differential equations. 

This problem was an attempt to study the performance of a particular fire control 
computer, replacing the human pilot by a suitable pllot-response-functlon. Thla same 
problem waa run on a REAC with a human pilot In the system. By observing the pilot's res¬ 
ponses over many REAC runs, Spangenberg and Banka were able to derive some sort of pilot- 
performance-function that should simulate the pilot's responses. Although the programming 
waa technically correct, the results were not what Spangenberg and Banka anticipated. 

More time was needed ln examining some of the logic used ln going from an analog-type 
simulation to a digital-type. However, since theses were due, further work could not be 
carried out. The conclusions that were drawn on thla Incomplete study can be found ln an 
MIT Flight Control Laboratory report |210-T20, "Response of the Human Pilot ln a Simulated 
Day-Superiority Type Fighter", by Spangenberg and Banka. 

C. Block 

Instrumentation 

Laboratory 

287 D. SAMPLED-DATA CONTACTOR SERVOMECHANISM 


283 B. A STUDY OF ERROR-REDUCTION IN INFORMATION SYSTEMS 


P.J.Arcldlacono 
Aeronautical Engineering 


- I W an attempt was maoe to simulate the error-reduction characterlatlea or a 15 -node 
network by the calculation of several arbitrary measures of redundancy Networks were 
choaen In which the measure of centrality aa outlined by Davelaa was constant. The mini¬ 
mum number of steps to completion of a task was also kept constant. Connections between 
nodes were represented digitally aa ones, no-conneotlona aa zeros. It was found 


The problem la to determine the response of a aumplea-data contactor servomech¬ 
anism which la disturbed by Gaussian nolne. The system considered has no Input other than 
the Gaussian nolae disturbance. Because of the noise, the system output at a given time la 
described by a list of possible system outputs with their associated relative probabilities. 
Thla Hat la oalled an output-distribution. Interest la centered on determining at sampling 
points the output-dlatrlbutton and/or lie variance. 

The output-dlatrlbutIon at sampling points la calculated by the following three 

general atepai 
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WHIRLWIND COPING A NO APPLICATION;' 


1) approximation)! are iada 10 reduce the number ol‘ possible system-outputs; 

2) the problem or calculating the output-distribution is stated In matrix form; 
the matrix proulem is solved on VNI. 

The three steps will ue discussed In order. 

T»ie number or poaoltile system-outputs la reduced by the following two methods: 

I) neglecting Improbable syate-.-outputs: 2) approximating the actual step response or the 
linear part or the system oy u step response wnlcn reaches steady-state orter a rew 
sampling periods. 

The output 1IstrlbutIon computation can ;e stated In matrix form In the following 
.runner. Consider a system-output vector V In a space whose coordinates are the possible 
System-outputs. The elements or V .ire then the probabilities or each or the possible system- 
output-". The trur.oformat Ion of the .• |. ,r ve toi or oi.e sampling point to Inc output vector 
or tic next sampling point can b< /level as a process or multiplying tne output vector bi¬ 
ll matrix of transition probabilities. For a syste- with no input other than a stationary 
Gaussian noise disturbance, the transition matrix Is constant. 

:a used to multiply the transition matri ■ - successive output vectors. 

The responses pt two abr.pleu-rtatu contactor ServomechBnlamB have been obtained on 
Kill. The Input to ! .1 output distribution vector, the tranaltlon 

matrix, snn the system-output corresponding to eacn vector, me procedure was coded by 
Q. Manoney of the MIT Digital Computer Uuoratory. Tne results or this problem have been 
Included In an S.M. thesis submitted by D. Chesler to the MIT Electrical Engineering Depart¬ 
ment . 

D. Chesler 

Electrical Engineering 

208 ATOMIC WAVE FUNCTION!: 

Electronic wave functions Tor atoms and uolecules (Tor states of arbitrary 
symmetry) -at. -e obtained to arbitrary accuracy o> -ethos wr.leh consists or three relat¬ 
ively Independent stages of calculation. 

Stage A - choice of some f 1. ;... set of a Ingle-electror. functions, a and evaluation 
of all possible one- and two-electron Integral; over tills set. The Integrals are matrix 
elements of operators occurring It. the any-electron Hamiltonian. 

Stage B - calculation or expansion coefficients or a set of orthonormal single- 
electron functions 

Stage C - calculation of configuration Interaction effects and resolution or 
degeneracies In the variational matrix for the Tuny-electron Hamiltonian. 


to calculate electronic wave functions, wltn extensive configuration Internet lor., for tin- 
r states or C, . 1 , and Or, 

H.K.Ieebet 
Solid State and 
Molecular Theory 
Group 

?*J C. HEAT TRANSFER THROUGH MIOll-SPEKD LAMINAR BOUNDARY 1AYERT 

In aeronaut lea, tne effects of tne vlacoalty and heat conductivity of air are 
generally conaldered to it confined -o u.ln layer of riulu, called the boundary laver. 
adjacent to the aircraft surface. tM* region, cerl-ili. terns In the Nuvlei-Stoker. 
equations and In tne en-rg .• equation for . viscous, neat-conductlng fluid may b- neglected 
•'id tne flow Is governed by tne , oundary la,er equations: 


<1> 

hi ,pu) * h lf,v| ’ 0 

(centlnulty) 

(2) 


(x-momeritum) 

(3) 

ft- 0 

(y-momentum) 

(A) 

"K*»8 *&*&<*& + *&* 

(energy) 


together with appropriate Boundary conditions at (he outer edge or the layer rind at me 
wall, an equation of state ana equations Wry* and fr as functions of the state ef the 
air. In these equations, applicable to two-dlmentlonal flows, the symbols are defined as 
follows: 

x co-ordinate In the W3ll along the direction of the free stream 

y co-ordinate normal to the wall 

u velocity component along x 

v velocity component along y 

p density 

p pressure 

h enthalpy 


If Stage a Is carried out, wltn certain me. tr lest Ions which are too complicated -o 
explain here. Stage C la greatly simplified. Techniques of group representation theory and 
perturoatlon metnoda cm be ISM. The actual .mount or calculation required Is verv small 
compared with Stages A and E. 

The essential difference between atomic mu molecular calculations Is that In the 
atomic case there exists a set of basic functions*| a , the analytic "later orbitals, which 
lend to rapid convergence In expansion or the self-conslstcnt jl orbitals, and for which all 
Integrals can he evaluated In closed form. ; 0 class of functions Is known which has both 
theoe propertl^o In tne molecular case. 


„ .. Pr 2? r " ' lrC •* all ‘|"} e Present which carry out stages A and b for atomic wave 

functions. These are not vet In rnetr most efficient form and some further programming Is 
noedeu to Join the. Into a single unit. Dr. '". K. Nesbet or the Solid State and Molecular 
Tr.roi". Group has progra. ei. t •• ralr.UaUon of Integrals, Stage A ror atoms, and the program 
lor transforming Integrals under iige h. These programs are uescrlueo In detail In the 
^quarterly Progress Report, olid state wr.u Molecular Theory Group, VIT, April, 195B. 

in this problem programs drveloped under Proolems lm., 2 10 and Z3 a are Mlng uscd 


Pr Prandtl numuer 

**■ vlacooity 

Exact solutions have seen obtained only for special classes of externa] flows 
and wall boundary conditions, isually with simplifying aasumptlons about the variation id' 
air properties. For practical problems Integral methodB are often used In which the 
equations are satisfied only on the average across tne boundary layer. The average 
values or -ass riux, momentum flux, energy flux, etc., are rrpresentrd h> ceri In i-o-i- 
ary laver thicknesses. Each thickness represents a ter In an equation o •.»Ine.l b: n. 
Integration of either equation (j) or (A) wit respect to y across the l«-.«r. followed by 
some manipulation. 
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" f * 

o 

■her* ( ), refers to conditions In the free stream and ( )„ to conditions at the wall. 

Theae quantities may be evaluated If the velocity and enthalpy profiles In the 
boundary layer and their relative scale, or thlcWneaa ratio, are given. It haa been round 
that the various profiles occurring In the Known exact solutions belong to one-parameter 
families - one Tor velocity and one Tor enthalpy - and It la generally assumed that the 
profiles in completely general flows belong to the same families. 

Four equations are needed In general to calculate the growth of a boundary layer 
with heat tranarer. Integral equations obtained from (5) and (4) describe the growth or 0 
and l If the velocity and enthalpy profiles are known. To obtain the two parameters 
defining the latter, two additional equations are needed. In some methods they are obtained 
from the original differential equations evaluated at the wall! In others new Integral equa¬ 
tions are used, obtained by first multiplying the differential equations by u or y and then 
integrating. 

An Integral method Is being developed at the Naval Supersonic laboratory which 
Involves essentially a new specification of the one-parameter ramlly of enthalpy profiles 
as a fixed profile for all flows with Isothermal walla on which a modifying profile may be 
superposed for riowa with non-Isothermal walla. The enthalpy profile parameter is then Just 
the ratio In which theae two profiles are combined. A new Integral equation had to be devel¬ 
oped for calculating the value of this parameter at any station along the surface from the 
given flow conditions. Actually two such equations were obtained, multiplying equation (4) 
first by u and then by y before Integrating, since It was not known in advance which would 
prove more suitable. Pour new thicknesses analogous to / e were Introduced: 

> 

o 


dy 


P(m,d) - J £ 

O 

T( "’ A> “ { s 7? 

u e®Pe 


h - h 


dy 


h - h 


^7 

h - h. 


) dy 


! i.'vC’ 4 ’ 


Olven either the Isothermal or the modifying enthalpy profile, each of these 
integrals (as well as J ) la a function of the velocity profile parameter, m, and the thick¬ 
ness ratio of the velocity and enthalpy profiles, 6. For each of 10 velocity profiles, val¬ 
ues were needed for 15 thickness ratios. A total of 1500 numerical Integrations were requir¬ 
ed. 


These Integrations were performed on the Whirlwind computer using the trapezoidal 
rule. The ten velocity profiles and two enthalpy profiles were fed Into the computer In 
tabular form. The tabular Interval of the velocity profiles was used aa the Interval for 
Integration. Since the thickness ratio varied, the value or the enthalpy In general had to 
be interpolated In the stored table. 

Tho Integration have been completed and work on the Integral method Is proceeding. 
The enthalpy profile modifying function was originally derived from the exact solution where 
the temperature potential across the boundary layer Is decreasing In the strearawlse direct¬ 
ion. This function may or may not be satisfactory for flows with Increasing temperature 
potential. If It la not, a new function will be derived, requiring a further 750 Integra¬ 
tions. 

J.A.P.H111 
Naval Supersonic 
Laboratory 
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290 N. POLARIZABILITY EPPECTS IN ATOMS AND MOLECULES 


The general problem la the calculation of the distortion produced on electronic 
wave functions for atoms under the Influence of a variety or applied electroatatlc fields 
and the subsequent use of this Information to construct aocurate electronic wove functions 
for molecules. There have been a number or programs written for Whirlwind which beer 
closely on this problem and the general aim is to utilize various combinations of these 
programs for production runs to obtain the large amount of numerical oats required. 

The work to date and that contemplated in the noar future has been solely con¬ 
cerned with the atomic part of the problem. For this problem the most easily manageable 
set of basis functions In which to expand the one-electron wave functions la 


r A+i e' ,r Yj (8.0) 

where A, a, m, f ere arbitrarily ohoaer. parameters. Various combinations or theae para¬ 
meters have been chosen and sets or basis functions have been built up to describe the 
polarization of the following atoms and Ions In a uniform field: 


H-. H e . L*, B«, B C *\ P-, N e . I,**, M g +Z , a/\ 


Cl", A, K*. end 0 a * Z . 

Similar seta of basis functions nave been set up to describe the distortion of the electro¬ 
nic wave function under the action of a point charge field for K*, Lj , and F". 

There are two stages In the calculation of these atomic wave functions and ener¬ 
gies which require extensive numerical work. There are first, the computation of atomic 
one- and two-electron integrals between the various basis functions and secondly, the 
algebraic transformation of these Integrals. The majority of the Integrals (those for the 
free atom) are being obtained from Dr. R. K. Nesbet's atomic Integral program which haa 
been described In these Progress Reports and in those of the Solid State and Molecular 
Theory Group. There are sufficiently few Integrals involving the perturbing field so that 
these have been done by hand computation and much of this work has been carried out by the 
R.L.E. Joint Computing Group. The second pert Involves a transformation program written 
by R.K.Nesbet (see Summary Report No. 41, p. 08, Problem 254) and the mechanization of the 
Roothaan scheme programmed by Dr. A. Meckler 6nd R.K.Nesbet (see description In previous 
reports). It is expected that runs using a combination of all these programs simultan¬ 
eously will soon be made. 


L.C.Allen 
Solid State and 
Molecular Theory 
Group 


291 B. DYNAMIC BUCKLING 

Thla problem la concerned with the response of structural columna to compressive 
end loadings which are applied very rapidly and which vary In time. The response may be 
elastic or plastic, the latter case being studied with the aid of VKI. The equations to 
be solved are a set of coupled ordinary differential equations of second order. The depend¬ 
ent variables are the coefficients of the normal modes of vibration of the structure, which 
are coupled In the plaetlc caae but not In the eleatlc oaae. The Kutta-Glll method Is being 
used to solve these equations. 

Two programs have been used so Tar, both being of a preliminary nature. They have 
served to acquaint the programmer with the machine and to define the spacing of the Inde¬ 
pendent variable which will give sufficiently accurate roaults. Tho second solution obtain¬ 
ed pointed out a problem of slow convergence which Is inherent in the mathematical work. 

The mathematics la presently being revised to correct thla trouble. 

R.E.Jones 
Civil Engineering 


292 A. COURSE 6.535, SPRIN0 1955 PRACTICE 


Course 6.555, "Introduction to Digital Computer Coding and 


c_a ...— 7fitc hi Mr HMn 


Logic" was given In the 
to solve one of the 
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following problems: 

1. program and oode a routine for the solution of a system of simultaneous 
linear equations by the Crout technique: 

2. code two routines for the solution of a differential equation by different 
numerical methods and compare the results; 

5. code a routine for the economization of a power series and compare the 

accuracy of the result with that or the power series trunoated to the same 
order. 

About 95 per cent of the students successfully completed the assignment. 


295 C. ROLLING BEARINGS 

The computer is to be used to solve the equations of motion of the rolling ele¬ 
ments. Eight transcendental equations Involving eight unknowns have been derived consid¬ 
ering applied thrust load, Inner ring speed, centrifugal force of rolling elements and 
normal and fractional forces In the contact zone. The direction and magnitude of the 
friction force In the contact zones are established by the relative slip, coefficient of 
friction and the local contact pressure. 

The program has been written and is in the testing stage. Extensive use Is made 
of library subroutines Including those for minimizing sum of squares by the gradient 
method, the Gauss quadrature method, and a new cube root routine. 

B. Rlakln 
Digital Computer 
laboratory 

294 C. WIND TUNNEL DATA REDUCTION 

The objective of the program was to see how tne Whirlwind computer oould be used 
to reduce data taken In the Naval Supersonic Wind Tunnel. In this particular example, 
forces and momenta (3 components) were measured cy strain gages mounted Inside the model. 
The strain gage output at present Is measured, In millivolts, by self balancing potentio¬ 
meters whose output In digital form Is automatically recorded on IBM punched cards. The 
data reduction process Involves converting the potentiometer readings to force and moment 
coefficients by multiplying by calibration constants and referring these coefficients to 
axea which are fixed in apace and do not move with the model. 

The program was oodod for the CS II computer with the output (run number, 
corrected angle of attack, dynamic pressure, base pressure coefficient, balance chamber 
pressure coefficient, drag coefficient, lift coefficient and pitching moment coefficient) 
recorded on magnetic tape. In addition, the machine was programmed to plot angle of attack, 
drag coefficient and pitching moment coefficient all as functions of lift coefficient. The 
desired results woro obtained for three runs. During each run, the angle of attack varied 
through about 10 values. 

Due to the method used In feeding data to the computer, unwanted tape numbers 
appeared between each row of tabulated data and on the plots so that a different method of 
storing the data would be used In future computations. Otherwise the program was success¬ 
fully completed and Indicated that the machine required about eight aeoonda per data point 
(a 'data point" being all of the data at one angle of attack). Although more efficient 
programing oould reduce this time slightly. It should still be a representative figure for 
more complicated programs (more components). All of the data from a good eight-hour shift 
of wind-tunnel time (about 20 runs, or 200 data points) could thus be reduced in about 
twenty-five minutes. 

L. Schlndel 
Naval Supersonic 
laboratory 


APPENDIX 

1. SYSTEMS ENGINEERING 


Revlalona In the Marginal Checking Procedures for the Drum s-. stem 

The programmed-marglna1-checklng racllltlea have been expanded to include addi¬ 
tional terminal equipment. Previously, this equipment was checked by taking itiat.ua 1 1-1111- 
glns during specially assigned maintenance periods. In order to Include u,e J.-m equip¬ 
ment In the dally programmed-marglnal-checklng routines, it was .eaasry to modify the 

marginal-checking equipment. A cheek program cycle for this type of terminal equipment 
waa necessarily long compared to the existing voltage variation cycle. It la now 
posslole to select the marginal-checking equipment It a npeolal mode writ-: will hold a 
preset excursion for an amount of time determined -y the program. 

A consolidated test program containing nine routines has been Written for the 
drum system. This program completely checks the equip art with a minimum amount or redun 
dance. Special care was taken In choosing program techniques and In assigning variation 
lines In an effort to reduce the amount of time necessary to perform the oheoka. The 
considerable reduction In time gained by this programming effort made It practical to 
Include the program as part of the dolly marginal-cheeking routine. The drum syatem, 
containing approximately ‘,500 cathodes, 1 b now checked In one half of the hour scheduled 
dolly for syatem maintenance. 

Analysis of Perlormance Records for the whirlwind Computer Syatem 

In September 195** the procedures for gathering and evaluating performance data 
on the computer syatem were somewhat revised. This was done to permit more comprehen¬ 
sive analyses of system reliability with particular emphasis on Interrupting failures. 

A description of these procedures as well an an analysis of the data gathered over a 20- 
week period ending 10 February 19b)) was published In Memorandum 6M-3410. 

Several riguroa are required to adequately describe the reliability of an elec¬ 
tronic system. In general, aystem reliability Is reflected In the amount or sched- led 
down time required for preventive maintenance. Since the amount of down time for differ¬ 
ent types of failures varies widely, the frequency of such failures Is also an important 
factor In deacrlblng aystom reliability. In the following table such reliability figures 
arc given. Results for the first 20-week period is reported In Memorsndum 6M-3410 are 
Hated along with similar data for the succeeding 5 months. 



20 Sept. 1954 
to 10 Feb. 1995 

11 Feb. 1955 
to 19 May 1955 

Total computer operating time, 

2675 hours 

1923 hours 

Total lost time 

92.7 hours 

65.2 hours 

Percentage operating time usable 

90.5 percent 

96.6 percent 

Average uninterrupted operating time 

Between failure Incidents 

10.6 hours 

10.9 hours 

Failure Incidents per 24-hour day 

2.15 

2.12 

Average lost time per Incident 

22.6 minutes 

23. minutes 

Average preventive maintenance time per day 

1.25 hours 

1.5 hours 
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2. PUBLICATIONS 


Project Whirlwind technical reporta and memoranoa are routinely dlatrlbuted 
to only a restricted group known to have a particular Interest In the Project, and to 
ASTIA (Armed Services Technical Information Agency) Document ervlee Center, Knott Build¬ 
ing, Dayton, Ohio. Requests for copies of Individual reports should Oe made to ASTIA. 


The following Is a Hat of memoranda published by the Scientific and Engineer¬ 
ing Computations Group during the past quarter. 


I.C.I., Ltd., England; Oulf Research and Development Co., Allegheny Refining Co.; 
Remington Rand Corp.; Meetinghouse Electric; Bell Telephone laboratories; The British 
Milk Marketing Board; The British Tabulating Co.; and Philips, N.V., Holland. 


No. 




DCL-66 

List of Memos on Programming and Coding 
for WWI 

A -A -55 

S and EC Oroup 

DCL-69 

S and EC Operator's Check List 

A-0-55 

J. Thompson 

DCL-71 

Notes on Tape Room Operating Procedure 

A-29-85 

M. Solomlta, 

H. Bello, and 

F. Shaw 

DCL-73 

Treatment of Plexo Coded Characters by The 
Comprehensive System (Revision of T^ble) 

8-5-55 

J. M. Prankovlch 

DCL-76 

The Availability of Floating Address Values 
Following CS Conversions 

5-20-55 

J. M. Prankovlch 

DCL-77 

Index of Subroutines Available for the 

WWI Computer 

6-1-55 

S and EC Oroup 


3. VISITORS 


Tours of the WWI Installation Include s showing or the film "Making Electrons 
Count", a computer demonstration, and an Informal discussion of the major computer com¬ 
ponents. During the past quarter, 9 groups totalling 188 people visited the computer 
Installation. Included In these groups were: 

April 28 

Newton High School 

May A 

M.I.T. Biology Department 

May 10 

Northeastern University 

May 11 

Operations Research Seminar 

May 26 

Class in “Introduction to Digital Computers - 
Coding and Logic", M.I.T. 

June 16 

Professor Gregory'a Summer oesslon class, M.I.T. 

June 2A 

Professor Hildebrand's Summer Session class, M.I.T. 


The procedure of holding Open House at the Digital Computer laboratory on the 
first Tuesday of each month has continued during this period. Three groups totalling 
58 persons visited the Laboratory at the Open House demonstrations. These persons re- 

E resented members and friends or the M.I.T. community, Boston Safe Deposit and Trust Co., 
lytheon, Oarlsnd Jr. College, Natick High School, laboratory for Insulation Research, 
John Hancock Life Insurance Co., Signal Manufacturing Co., and the Massachusetts Memorial 
Hospitals. 

During the past quarter there were also 78 Individuals who made brief tours of 
the computer Installation at different times. Some of the companies represented by 
these Individuals are Chanee-Vought Aircraft; Bridgestone Tire Co. and Oktsu Tire Co. of 
Japan; Godfrey L. Cabot Co.; Institute de Statlstlque, Paris; British Grain and Seed Co.- 
U. S, Army Signal Corps; Whirlpool Co., Link Aviation Co.; Monsanto Chemical Co.; 
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